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SECTION  1 
INTRODUCTION 

Current  research  efforts  in  ballistic  reentry  vehicle  aero¬ 
dynamics  are  prim?rily  concerned  with  the;  improvement  of  vehicle  targeting 
accuracy.  Accurate  evaluation  of  possible  targeting  errors  requires  a 
detailed  understanding  of  all  mechanisms  that  may  deflect  the  vehicle 
away  from  its  nominal  ballistic  trajectory.  Of  the  dispersion  errors 
that  can  be  attributed  directly  to  the  reentry  vehicle,  the  low  altitude 
roll -trim  effect  is  one  of  the  prime  contributors  to  miss  distance. 

Roll -trim  dispersion  results  when  normal  forces,  such  as 
created  by  a  trim  angle  of  attack  condition,  are  not  integrated  cut  by 
the  spin  of  the  vehicle.  The  characterization  of  such  dispersion  re¬ 
quires  a  coupling  of  the  vehicle's  dynamics  with  its  aerodynamic  charac¬ 
teristics  along  the  entire  entry  trajectory.  This  effort  is  aimed  at  ex¬ 
tending  current  aerodynamic  prediction  capabilities  relative  to  the  low 
altitude  roll -trim  dispersion  problem. 

The  current  generation  of  ballistic  reentry  vehicles  are 
typically  slender  blunted  cones  or  biconic  configurations,  with  the  nose- 
tips  generally  being  fabricated  from  woven  carbonaceous  materials.  During 
reentry  the  severity  of  the  aerothermodynamic  environment  causes  ablation 
of  the  nosetip  material,  leading  to  both  axial  recession  of  the  nosetip 
and  to  alterations  in  the  basic  shape  of  the  nose. 

At  higher  altitudes  (100  KFT  i  h  i  50  KFT)  the  flow  in  the 
nosetip  boundary  layer  remains  laminar,  resulting  primarily  in  blunting 
of  the  nosetip.  Below  approximately  50  KFT,  however,  as  the  nosetip 
boundary  layer  is  passing  through  transition  to  a  fully  turbulent 
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state,  the  increased  heating  levels  lead  to  a  sharpening  of  the  nosetip 
shape,  as  illustrated  in  Figure  1.1. 

Because  of  circumferential  variations  in  the  onset  and  pro¬ 
gression  of  nosetip  transition,  asymmetric  nose  geometries  can  result. 
The  mechanisms  governing  this  transition  process,  such  as  surface  rough¬ 
ness  variations, are  generally  evaluated  statistically,  as  by  Dirling1  . 


Figure  1.1.  Typical  Nose  Shape  Progression 
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The  development  of  an  asytrmetric  nosetip  shape  on  an  otherwise 
axi symmetric  body  will  lead  to  the  development  of  a  trim  angle  of  attack 
and  corresponding  side  forces  which  tend  to  deflect  the  vehicle  from  its 
ballistic  trajectory.  To  minimize  these  trim  effects,  reentry  vehicles 
are  spun  prior  to  their  reentering  the  atmosphere;  thus  the  integrated 
resultant  of  a  nody-fixed  lift  force  over  one  revolution  of  the  body  will 
be  nearly  zero.  However,  rapid  variations  in  nose  shape  and  angle  of 
attack  (and  hence  lift  force)  and  roll  rate  can  result  in  a  non-zero 
resultant  force,  leading  to  trajectory  deflection  due  to  roll-trim  dis¬ 
persion. 

Because  of  the  inherent  uncertainties  in  nose  shape  change 
predictions,  the  roll-trim  effect  is  usually  evaluated  statistically, 
as  by  Pettus,  Larmour,  and  Palmer  2  .  Given  a  nose  shape,  however,  the 
evaluation  of  aerodynamic  characteristics,  a  necessary  part  of  any  roll- 
trim  evaluation,  is  a  deterministic  problem. 

Aside  from  expensive  and  time-consuming  wind  tunnel  tests,  the 
most  accurate  and  reliable  method  for  the  prediction  of  aerodynamic 
characteristics  is  the  numerical  integration  of  the  invsicid  equations 
of  fluid  motion.  Fortunately,  at  the  flight  conditions  of  interest  (and 
in  most  hypersonic  wind  tunnels  simulating  reentry  conditions  at  low 
altitudes)  the  Reynolds  number  is  sufficiently  large  that  the  shock 
layer  is  almost  entirely  inviscid,  except  for  the  thin  boundary  layer 
adjacent  to  the  vehicle  surface.  Additionally,  the  flow  is  in  the  weak 
interaction  regime,  as  defined  in  Hayes  and  Probstein3  ,  where  viscous 
shear  and  induced  pressure  effects  significantly  affect  only  the  axial 
force  experienced  by  the  vehicle.  Other  vehicle  forces  and  moments 
(normal  and  side  forces,  pitching  and  yawing  moments)  can  then  be 
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accurately  determined  solely  through  consideration  of  the  inviscid 
pressure  distribution. 

In  the  past  decade  many  numerical  procedures  have  been  developed 
for  the  calculation  of  inviscid  aerodynamic  characteristics  for  ballistic 
reentry  vehicles.  These  techniques  have  proven  to  be  valuable  adjuncts 
to  the  design  process  and  have,  to  some  extent,  lessened  the  ne.:.d  of 
performing  extensive  wind  tunnel  tests  to  validate  proposed  configurations. 
The  three-dimensional  numerical  procedures  currently  in  use  consist  of 
two  parts:  a  transonic  flow  field  procedure  to  treat  the  subsonic  region 
surrounding  tne  stagnation  point,  and  an  afterbody  procedure  to  treat  the 
the  supersonic  flow  on  the  vehicle  frustum. 

The  inviscid  afterbody  flow  field  problem  is  now  well  in  hand 
for  the  simple  axisymmetric  frusta  found  on  ballistic  reentry  vehicles. 

(In  addition,  ballistic  vehicles  at  low  altitudes  generally  do  not 
develop  large  angles  of  attack  which  would  lead  to  flow  separation  on 
the  leeside  of  the  afterbody,  invalidating  the  inviscid  assumption.)  The 
existing  inviscid  transonic  nosetip  flow  field  capability  is  restricted, 
however,  to  convex  shapes,  where  strong  embedded  shock  waves,  such  as 
shown  in  the  Schlieren  photograph  in  Figure  1.2,  do  not  occur.  Further¬ 
more,  other  restrictions  arise  even  for  convex  shapes,  when  the  coordinate 
system  used  in  the  calculation  is  not  closely  aligned  with  the  surface  of 
the  nosetip.  (These  shortcomings  of  the  current  techniques  were  identi¬ 
fied  by  Hall,  Kyriss,  Truncellito,  and  Martel! ucci  4  .) 


12. 


of  Ablated  Nosetip  with  Embedded  Shock 


The  goal  of  the  current  research  effort  is  to  eliminate  the 
above  two  restrictions  on  inviscid  transonic  flow  field  techniques.  By 
extending  the  range  of  nosetip  shapes  that  can  be  analyzed  numerically 
to  include  slender  and  indented  shapes  such  as  have  been  observed  in 
flight,  the  capability  for  accurate  evaluation  of  roll-trim  dispersion 
will  be  greatly  expanded.  In  addition  this  nrw  capability  will  allow 
more  accurate  nose  shape  reconstruction  efforts  (in  which  a  nose  shape 
is  sought  that  produces  aerodynamic  characteristics  that  agree  with 
those  derived  from  the  vehicle  motion  observed  in  flight),  as  described 
by  Hall  and  Nowlan5  .  Furthermore,  this  transonic  flow  field  technique 
will  be  applicable  to  maneuvering  as  well  as  ballistic  reentry  vehicles, 
since  autopilot  design  for  maneuvering  vehicles  must  account  for  the 
aerodynamic  characteristics  that  result  from  ablated  nosetip  geometries. 

The  approach  taken  to  eliminate  these  deficiencies  of  the 
current  nosetip  flow  field  procedures  is  outlined  in  Section  2.2. 

Section  3.0  details  the  conformal  mapping  transformation  used  to  produce 
a  coordinate  system  closely  aligned  with  the  body  surface,  and  Section 
4.0  describes  the  procedures  used  for  the  calculation  of  embedded  shocks. 
In  Section  5.0,  details  are  provided  on  the  numerical  procedures  used  in 
this  new  transonic  flow  field  technique,  which  is  validated  by  compari¬ 
sons  to  data  in  Section  6.0. 


14. 


•  / 


SECTION  2 


PROBLEM  DEFINITION 

2.1  STATEMENT  OF  THE  PROBLEM 

The  problem  being  examined  in  this  research  effort  involves 
the  numerical  prediction  of  the  inviscid  aerodynamic  characteristics 
o *  ballistic  reentry  vehicles  with  asymmetric,  ablated  nosetips.  In 
particular,  emphasis  is  placed  on  the  development  of  a  numerical  tech¬ 
nique  to  determine  the  inviscid  flow  field  about  a  three-dimensional, 
asynmetric  nosetip  1r  a  uniform  hypersonic  or  supersonic  freestream  flow. 
The  determination  of  the  nosetip  flow  field  is  a  necessary  first  step 
for  the  prediction  of  total  vehicle  aerodynamics. 

The  assumption  is  made  in  this  analysis  that  inviscid  flow 
theory  is  adequate  to  accurately  predict  the  aerodynamic  characteristics 
of  reentry  vehicles  at  altitudes  where  asynwstric  nose  shapes  can  result 
from  ablation  (h  s  50  KFT).  (Accurate  calculation  of  drag  forces  will 
also  require  consideration  of  surface  shear  and  induced  pressure  effects, 
however.)  Implicit  in  the  assumption  of  inviscid  flow  are  the  require¬ 
ments  that  the  thin  boundary  layer  assumption  be  valid  and  that  no  regions 
of  separated  flow  exist  on  the  vehicle. 

Moretti  and  Salas  6  ,  in  their  analysis  of  viscous  rarefied 
flows,  have  presented  a  breakdown  of  the  various  flow  regimes  that  might 
be  expected  as  a  function  of  freestream  Mach  and  Reynolds  numbers  (with 
the  Reynolds  number  based  on  nosetip  radius  for  a  spherical  nose),  as 
depicted  in  Figure  2.1.  Also  indicated  in  this  figure  is  an  M^  -  Re^ 
history  for  a  typical  modern  ballistic  reentry  vehicle  as  a  function  of 
altitude.  Defining  the  thin  boundary  layer  regime  as  the  region  where 
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REYNOLDS  NUMBER  Re^  = 

Reoo-Mco  Map  of  Flow  Regimes  with  Typical  Reentry  Vehicle 
Trajectory  Superimposed 


the  boundary  layer  thickness  (5)  Is  less  than  1%  of  the  shock  layer 
thickness  (d),  this  figure  clearly  indicates  the  validity  of  the  in- 
viscid  assumption  for  the  problem  being  considered  here. 

The  hypersonic  flow  over  a  blunt  nosetip  is  characterized  by  a 
detached  bow  shock  wave  separating  the  shock  layer  from  the  undisturbed 
freestream  flow.  In  the  vicinity  of  the  stagnation  point,  where  the 
bow  shock  is  nearly  normal  to  the  freestream  velocity  vector,  the  flow 

in  the  shock  layer  is  locally  subsonic;  thus,  the  steady  flow  problem  in 

\ 

this  region  has  an  elliptic  character.  As  the  bow  shock  curves  back 
around  the  body  and  becomes  more  oblique  to  the  freestream  velocity 
vector,  and  as  the  shock  layer  flow  expands  around  the  nose,  the  flow 
becomes  locally  supersonic,  and  the  steady  flow  problem  becomes  hyper¬ 
bolic.  Other  complications  can  arise  in  this  basic  inviscid  flow  field 
structure  if  the  body  surface  has  indented  regions  producing  embedded 
shock  waves.  Depending  on  the  shock  strength,  the  flow  behind  such  an 
embedded  shock  could  be  either  subsonic  or  supersonic. 

Because  of  this  variety  of  flow  conditions  that  can  be  en¬ 
countered  in  the  blunt  body  problem,  it  is  convenient  to  seek  the  steady 
solution  as  the  asymptotic  limit  of  the  time-dependent  problem,  since 
the  unsteady  flow  equations  are  hyperbolic  in  time,  regardless  of  the 
local  Mach  number.  Furthermore,  since  the  location  of  the  bow  shock 
wave  is  unknown  a  priori  and  must  be  determined  as  part  of  the  solution 
procedure,  this  time -asymptotic  approach  has  the  additional  benefit  of 
allowing  the  calculation  of  the  time  history  of  the  shock  shape  starting 
from  an  assumed  initial  shock  position. 

The  numerical  procedure  to  be  used  in  the  solution  of  this 
time-dependent  problem  is  an  explicit,  second-order  accurate  finite- 


difference  p^jceoure.  Similar  schemes  have  been  developed  previously, 
and  are  cuirently  in  wide  use;  however,  these  procedures  are  limited  in 
their  ability  to  treat  slender  ablated  nosetlp  shapes,  and  in  their 
ability  to  treat  embedded  shocks.  This  research  effort  is  devoted  to 
the  development  of  a  procedure  that  eliminates  the  deficiencies  observed 
in  other  transonic  time-dependent  codes.  In  particular,  this  requires 
development  of  a  generalized  coordinate  system  that  is  capable  of  being 
closely  aligned  with  the  three-dimensional  body  surface  for  abritrary 
body  geometries,  as  well  as  the  development  of  a  procedure  for  the  cal¬ 
culation  of  three-dimensional  embedded  shocks  on  indented  nosetip  shapes. 

Coupling  this  nosetip  transonic  procedure  to  an  existing 
supersonic  afterbody  code  will  thus  allow  accurate  theoretical  assess¬ 
ment  of  the  effect  of  ablated  nosetip  geometries  on  the  performance  of 
the  total  vehicle  for  many  nose  shapes  that  could  net  previously  be 
analyzed. 

2.2  OUTLINF  OF  APPROACH 

The  fundamental  approach  selected  in  this  research  effort  for 
the  solution  of  the  blunt  body  problem  for  ablated  asymmetric  nosetips 
is  the  time -dependent  relaxation  approach.  This  technique  has  been  widely 
used  in  previous  work  and  has  several  advantages  directly  related  to 
difficulties  associated  with  the  blunt  body  problem.  In  particular,  the 
time-dependent  approach  allows  the  use  of  a  convenient  forward-marching 
(in  time)  numerical  algorithm,  avoiding  many  of  the  difficulties  other¬ 
wise  encountered  in  the  steady  flow  problem. 

The  numerical  scheme  selected  for  this  analysis  is  in  many  ways 
similar  to  that  used  in  other  procedures.  For  example,  the  treatment  of 
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field  points  In  this  algorithm  (when  no  embedded  shocks  are  present)  Is 
based  on  the  second-order  accurate  explicit  MacCormack7  predictor- 
corrector  finite  differencing  scheme.  This  particular  scheme  has  found 
wide  application  to  computational  fluid  dynamics  problems  (e.g.,  References 
8,  9,  10,  11,  12)  because  of  Its  high  degree  of  accuracy  and  ease  of 
Implementation.  Boundary  points  at  the  body  surface  and  at  the  outer 
(bow)  shock  are  treated  using  the  Ken tzer -Moretti  predictor-corrector 
boundary  point  procedure,  in  which  a  discretization  of  boundary  conditions 
suggested  by  Kentzer13  was  extended  to  a  predictor-corrector  format  by 
Moretti  and  Pandolfi9.  This  boundary  point  procedure  has  found  wide 
application  in  computational  fluid  dynamics,  and  its  properties  have  been 
analyzed  and  discussed  by  Hall14. 

A  conformal  mapping  technique  was  selected  to  define  the  coordi¬ 
nate  system  for  this  problem  because  of  its  ability  to  preserve  local 
angles  under  the  transformation.  Thus,  by  formulating  a  mapping  in  which 
the  image  of  the  body  surface  is  a  nearly  horizontal  line,  and  selecting 
another  coordinate  direction  to  be  the  vertical  direction  in  the  trans¬ 
formed  space  (and  hence  nearly  normal  to  the  body  image),  the  resulting 
grid  in  physical  space  will  consist  of  surfaces  closely  aligned  with  and 
nearly  normal  to  the  body  surface.  The  ability  to  automatically  generate 
such  a  coordinate  system  for  ablated  asymmetric  nosetip  shapes  is  critical 
to  the  success  of  the  numerical  algorithm  in  computing  inviscid  flows 
about  such  shapes. 

The  coordinate  system  developed  in  this  research  is  based  on 
the  "hinge  point"  concept  of  Moretti,  as  developed  in  References  15,  16, 
and  17.  Application  of  this  technique  to  the  asymmetric  nosetip  problem 


19. 


has  required  the  extension  of  this  technique  to  three  dimensions;  this 
development  is  described  In  Section  3.2. 

To  treat  the  embedded  shock  problem  in  this  research  effort, 
the  shock-capturing  approach  has  been  selected,  in  which  the  structure  of 
the  embedded  shock  is  approximated,  but  for  which  no  special  logic  is 
required  to  explicitly  treat  the  shock.  Two  shock-capturing  approaches 
are  examined:  the  conservation  formulation,  discussed  in  Section  4.1, 
and  the  X-differencing  approach,  discussed  in  Section  4.2.  Axisymmetric 
versions  of  both  of  these  techniques  have  been  developed,  and  a  comparison 
of  the  results  of  these  two  approaches  is  made  in  Section  4.3.  Based  on 
these  comparisons,  it  is  concluded  that  the  X-differencing  approach  is 
the  superior  method,  and,  accordingly,  is  extended  to  three-dimensions  in 
Section  4.4. 

Axi symmetric,  inviscid,  time-dependent  procedures  with  a  shock¬ 
capturing  approach  to  the  treatment  of  embedded  shocks  have,  of  course, 
been  developed  previously,  notably  by  Kutler,  Chakravarthy,  and  Lombard18 
(using  the  conservation  form)  and  by  Moretti18  (using  the  X-differencing 
approach).  The  successful  application  of  the  X-shock-capturing  technique 
to  the  three-dimensional  time-dependent  embedded  shock  problem  is,  however, 
new. 

The  final  portion  of  this  effort  is  devoted  to  the  validation 
of  the  resulting  numerical  technique  for  the  calculation  of  inviscid  nose- 
tip  flow  fields.  For  simple  nosetip  shapes  (e.g.,  spheres)  this  new 
technique  is  compared  to  proven  flow  field  codes,  such  as  that  developed 
by  Kyriss  and  Harris®.  For  other  shapes,  representative  of  nosetip  shapes 
that  result  from  the  ablation  process,  comparisons  are  made  with  wind 
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tunnel  measurements  of  surface  pressure  and  bow  shock  shape.  Un¬ 
fortunately,  most  of  the  existing  wind  tunnel  data  providing  these  details 
on  nosetip  flows  are  available  only  for  axi symmetric  shapes,  as  in 
Reeves,  Todlsco,  Lin,  and  Pall  one2 9  and  Jackson  and  Baker21. 

A  large  body  of  wind  tunnel  data  does  exist  for  the  total  aero¬ 
dynamics  of  reentry  vehicles  with  asymmetric  nosetips.  Coupling  the  new 
nosetip  flow  field  procedure  to  an  existing  afterbody  code,  comparisons 
are  made  between  predicted  and  measured  vehicle  forces  and  moments,  thus 
providing  an  indirect  means  of  verifying  the  accuracy  of  the  nosetip 
calculation.  These  comparisons  are  presented  in  Section  6.4. 
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SECTION  3 

COORDINATE  SYSTEM  AND  GOVERNING  EQUATIONS 


This  section  provides  details  on  the  coordinate  transformation 
developed  in  this  effort  that  is  capable  of  mapping  the  surface  ot 
ablated,  asymmetric  nosetip  geometries  onto  a  nearly  horizontal  surface, 
thus  producing  a  coordinate  grid  cl)s?ly  aligned  with  the  body  geometry. 
This  mapping  is  then  used  to  generate  the  three-dimensional  inviscid 
time-dependent  equations  of  fluid  motion  written  in  terms  of  the  new 
coordinates.  A  final  computational  transformation  is  described  that 
maps  the  transformed  shock  layer  onto  a  regular,  equally  spaced  grid. 

The  equations  derived  in  this  section  are  written  in  non¬ 
conservation  form;  i.e.,  the  dependent  variables  are  the  primary  flow 
variables.  This  form  of  the  governing  equations  is  the  form  used  for 
flow  calculations  when  embedded  shocks  are  not  present,  and  as  the  basis 
of  the  X-differencing  shock-capturing  scheme  described  in  Section  4.2. 

The  conservation  form  of  the  governing  equations  is  discussed  in  Section 
4.1. 

3.1  INVISCID  EQUATIONS  OF  MOTION 

The  three-dimensional  time-dependent  inviscid  equations  of 
motion  in  non-conservation  form  may  be  written  in  cylindrical  coordinates 
as 

Pt  +  UPX  +  VPy  +  WP^/y  +  y(Ux  +  Vy  +  W^/y  +  V/y)  =  0  (3.1) 

Ut  +  UUX  +  VUy  +  WU^/y  +  pPx/p  =  0  (3.2) 

Vt  +  UVX  +  VVy  +  WV^/y  -  w2/y  +  pPy/P  =  0  (3.3) 
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(3.4) 


Wt  +  UWX  +  VW  +  WW^/y  +  VW/y  +  pP^/py  =  0 
st  +  Usx  +  Vsy  +  Ws<jA  “  °  (3.5) 

where 

$  =  UI  +  VJ  +  WK 

A  A  /V 

with  I,  J,  and  K  being  the  unit  vectors  in  the  x,  y,  and  <j>  directions, 
respectively.  (In  this  cylindrical  system,  x  is  the  axial,  y  the  radial, 
and  <J>  the  circumferential  coordinate.)  In  this  formulation  the  dependent 
thermodynamic  variables  are  P  and  s,  where 

P  =  in  p  (3.6) 

and  s  is  some  suitable  analog  of  the  entropy.  The  choice  of  P  as  a 
dependent  variable  is  motivated  by  computational  considerations,  since 
the  logarithm  of  pressure  throughout  the  shock  layer  will  not  vary  by 
several  orders  of  magnitude  as  the  pressure  might;  thus,  one  can  expect 
more  accurate  finite  difference  representations  of  derivatives  of  P  than 
could  be  expected  for  p. 

For  closure  of  this  system  of  partial  differential  equations, 
a  thermodynamic  equation  of  state  of  the  form 

P  =  p(p,s)  (3.7) 

is  required.  For  an  equilibrium  real  gas  calculation,  the  relation  em¬ 
bodied  in  Equation  (3.7)  may  be  provided  either  through  tabulations  of 
the  thermodynamic  properties  or  through  an  appropriate  curve  fit  of  the 
thermodynamic  data.  In  the  case  of  a  thermally  and  calorically  perfect 
(ideal)  gas,  this  thermodynamic  relation  may  be  expressed  implicitly  as 


S  ■  (in  p  -Y  *npy(Y-1) 


(3.8) 


where  the  thermodynamic  variable  s  is  defined  in  terms  of  the  entropy 
?  as 

s  =  (?  -  sJ/R  (3.9) 

with  y  the  isentropic  exponent  and  R  the  gas  constant.  Inversion  of 
Equation  (3.8)  yields 

Ills 

p  =  pl/Ye  Y  .  (3.10) 

Finally,  to  complete  the  definition  of  the  mathematical  problem, 
initial  and  boundary  conditions  must  be  specified.  Since  the  steady-state 
solution  is  sought  as  the  asymptotic  limit  of  the  transient  problem,  the 
specification  of  an  initial  flow  field  is  required.  Details  on  the  de¬ 
finition  of  this  assumed  initial  flow  field  are  presented  in  Section  5.1. 

Boundary  conditions  for  this  problem  must  be  specified  at 
the  boundaries  of  the  region  being  computed:  at  the  bow  shock  wave 
y  =  ys(x,<|>),  on  the  body  surface  y  ■  y^(x ,<|)) ,  and  on  some  downstream 
boundary,  running  between  the  body  and  the  shock.  The  location  of  the 
downstream  boundary  is  arbitrary,  subject  only  to  the  restriction  that 
the  flow  across  this  boundary  be  supersonic.  As  long  as  this  boundary 
is  entirely  supersonic,  no  condition  need  be  imposed  there,  since  the 
range  of  influence  of  this  boundary  will  then  not  extend  back  into  the 
region  being  computed. 

At  the  bow  shock,  whose  position  is  unknown  a  priori  and 
must  be  determined  as  part  of  the  solution  procedure,  the  appropriate 
boundary  conditions  are  given  by  the  familiar  Rankine-Hugoniot  conditions. 


By  incorporating  differential  forms  of  these  relations  for  conservation 
of  mass,  momentum,  and  energy  across  the  shock  into  a  characteristic 
compatibility  relation,  an  equation  for  the  shock  acceleration  is  obtained, 
which  may  be  integrated  to  yield  shock  velocity  and  position.  This 
numerical  scheme  is  described  in  more  detail  in  Section  5.4. 

At  the  body  the  appropriate  boundary  condition  to  be  imposed 
is  the  inviscid  kinematic  boundary  condition,  which  requires  that  there 
be  no  velocity  component  normal  to  the  body  surface.  This  condition  is 
applied  in  conjunction  with  a  characteristic  compatibility  condition  to 
develop  a  numerical  procedure  for  body  points  as  described  in  Section  5.3. 

Also  of  interest  in  the  treatment  of  boundary  conditions  for 
this  problem  is  the  value  of  entropy  that  applies  along  the  streamline 
that  wets  the  body  surface.  It  is  frequently  assumed  that  the  surface 
entropy  for  inviscid  flows  is  exactly  the  normal  shock  value  of  entropy, 
but  this  can  be  proven  only  for  axi symmetric  flows.  Numerical  results 
of  Swigart22  and  the  experimental  results  of  Xerikos  and  Anderson23 
indicate  that  this  assumption  may  not  be  true  and  that  the  normal  shock 
streamline  does  not  wet  the  body  surface.  However,  in  his  survey  paper, 
Rusanov24  argues  that  the  results  of  Swigart's  calculations  using  an 
inverse  method  are  inconclusive  because  of  the  inherent  assumptions  and 
computational  errors.  Additionally,  Rusanov  points  out  that  his  own 
computational  results  using  a  finite  difference  procedure  produced 
variations  between  the  computed  surface  entropy  and  the  normal  shock 
entropy  of  less  than  0.1%,  which  is  within  the  error  level  of  his  calcu¬ 
lation.  From  the  examination  of  his  studies  and  the  results  of  others, 
Rusanov  concludes  that  there  is  no  firm  evidence  of  the  surface  entropy's 


not  having  the  normal  shock  value,  although  there  is  likewise  no 
proof  that  these  values  coincide. 

From  a  practical  standpoint,  the  question  of  the  value  of  the 
surface  entropy  is  not  critical,  since  even  the  variations  in  surface 
entropy  claimed  by  Swigart  and  Xerikos  and  Anderson  produce  only  small 
perturbations  on  the  other  flow  variables  (e.g.,  density,  velocity). 
Accordingly,  the  body  surface  is  assigned  the  known  normal  shock  value 
of  entropy  in  this  problem. 

Circumferentially,  the  boundary  condition  to  be  imposed  in 
this  problem  is  that  of  periodicity;  i.e.,  the  solution  at  <(>  =  0  must 
coincide  with  the  solution  at  <j>  -  2ir.  For  the  case  of  a  pitch  plane  of 
symmetry  (geometric  symmetry  and  no  sideslip),  the  calculation  need  be 
performed  only  from  d>  *  0  to  $  *  it,  and  the  circumferential  boundary 
conditions  simply  require  synmetry  about  the  pitch  plane. 

3.2  THREE-DIMENSIONAL  CONFORMAL  TRANSFORMATION 

A  major  portion  of  this  research  is  devoted  to  the  development 
of  a  generalized,  three-dimensional  coordinate  transformation  that  is 
capable  of  producing  a  coordinate  surface  closely  aligned  with  the  body 
surface.  As  outlined  earlier,  the  method  that  has  resulted  is  based  on 
an  idea  of  Moretti's17  for  axisymmetric  time-dependent  calculations. 

The  general  coordinate  transformation  used  takes  the  functional 

form 


£  ■  £(x,y,<t>) 

(3.11) 

n  =  n(x,y,<!>) 

(3.12) 

d  =  $ 

(3.13) 

T  =  t 

(3.14) 
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which  Implies  that  the  transformation  of  the  spatial  coordinates  is 
independent  of  time.  Furthermore,  <J>  *  constant  planes  are  transformed 
directly  to  0  *  constant  planes,  thus  retaining  a  somewhat  "cylindrical" 
quality  to  the  transformation. 

Prior  to  reentry,  ballistic  vehicles  are  initially  axisym- 
metric,  and  It  may  thus  be  expected  that  ablated  asymmetric  nosetip 
shapes  that  develop  during  reentry  will  retain  some  "axisymr' trie" 
character.  In  other  words,  since  the  <J>  *  constant  planes  wil  be  normal 
to  the  vehicle  surface  prior  to  reentry,  it  is  reasonable  to  expect  that 
the  simple  transformation  9  *  <p  will  lead  to  8  =  constant  planes  that  are 
nearly  normal  to  the  surface  of  the  ablated  nosetip,  even  though  the 
ablated  shape  may  not  be  truly  axi symmetric. 

Within  each  <(>  =  constant  plane,  then,  the  transformation  re¬ 
duces  to  the  form 

5  =  ?(x,y)  (3.15) 

n  ■  n(x,y)  .  (3.16) 

Since  it  is  desirable  to  have  a  coordinate  grid  closely  aligned  with 
the  body  geometry  (and  hence  with  the  streamlines  of  the  flow),  a 
transformation  is  sought  that  closely  aligns  the  5  direction  with  the 
body  surface  (within  a  <j>  =  constant  plane).  In  order  to  have  the 
n-.di  recti  on  normal  to  the  5-di  recti  on  at  all  points  (and  hence  nearly 
normal  to  the  body  surface),  a  conformal  transformation  is  sought,  since 
under  a  conformal  transformation,  the  orthogonal  (x,y)  grid  maps  onto 
an  orthogonal  (5,n)  grid. 


27. 


Conformal  transformations  from  the  *  x  +  1y  space  to 
the  c  ■  £  +  in  space  can  then  be  developed  independently  in  each  $-plane. 
These  transformations  rely  on  the  concept  of  "hinge  points"  as  developed 
by  Moretti  15,16,17  to  ensure  that  the  ^-direction  is  closely  aligned  with 
the  body  surface. 

The  concept  behind  this  "hinge  point"  approach  is  to  define 
a  sequence  of  points  in  the  z-j  space  that  lie  close  to  the  body  surface 
and  define  an  approximate  equivalent  body  shape.  A  sequence  of  con¬ 
formal  transformations  is  then  applied  to  map  each  of  these  hinge  points 
in  turn  onto  the  horizontal  axis;  if  the  hinge  points  in  the  z-j  space 
accurately  simulate  the  body  contour,  the  resulting  transformed  contour 
will  then  be  nearly  horizontal  (i.e.,  will  be  closely  aligned  with  a 
coordinate  surface). 

for  the  mapping  function  developed  in  this  research,  which 
has  been  adopted  by  Moretti17  for  axi symmetric  calculations,  hinge  points 

i.L 

are  defined  as  illustrated  in  Figure  3.1.  Let  h.  .  .  denote  the  i 

■  *  J  *  K 

hinge  point  in  the  transformed  space  (j  =  1  is  (x,y)  space)  in  the 
4)  =  <t>k  plane.  It  is  required  that  h^  k  be  located  on  the  nosetip 
centerline  outside  of  the  body  and  that  h„  ,  t  be  located  on  the  center- 
line  inside  the  body.  The  remaining  hinge  points  h.  1  k,  i  *  3,4,. ...JC 
are  selected  so  as  to  model  the  body  contour.  Note  from  Figure  3.1  that 
this  specification  of  points  produces  OA  =  JC  -  2  "corners"  which  must 
be  eliminated  by  the  mapping  sequence  to  have  all  hinge  point  images  on 
the  horizontal  axis  (in  the  transformed  space). 


To  eliminate  each  corner  in  succession,  the  mappings,  which 
have  been  developed  as  part  of  this  effort,  of  the  form 


zj+l,k  '  1  s  tzj,k  "  Vl.J.k1 


5*  i. 

J,k 


(3.17) 


are  applied  sequentially  for  j  =  1,  2,...,JA.  The  form  of  this  transforma¬ 
tion  is  related  to  the  Schwarz-Christoffel  transformation  and  indeed  may 
be  regarded  as  a  "point-wise"  Schwarz-Christoffel  transformation.  By 
proper  selection  of  the  exponents  6.  ,,,  defined  from 

J  » K 


6j,k 


rr-tan-1 

Rethj+2,j,k  "  Vl,i,k] 


(3.18) 


these  mappings  have  the  required  property  of  maintaining  all  hinge 
points  h^j^,  i  s  j  on  the  real  axis,  while  mapping  h^  ^  k  onto  the 
real  axis. 


The  application  of  this  mapping  is  illustrated  in  Figure  3.2, 
showing  how  each  of  the  OA  corners  is  eliminated  successively,  resulting 
in  all  h..  jg  k  (with  JB  =  JA  +  1  =  JC  -  1)  lying  on  the  real  axis  in 
the  ZjB>k  space.  It  is  important  to  note  that  straight  line  segments 
between  hinge  points  in  the  z^  k  space  are  not  maintained  as  straight 
segments  under  this  sequence  of  transformations.  Since  each  of  these 
intermediate  transformations  is  conformal,  the  sequence  of  mappings  will 
itself  be  a  informal  transformation. 
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Two  further  transformations  are  required  to  complete  the 
specification  of  a  suitable  coordinate  grid.  First,  it  Is  beneficial 
to  have  the  transformed  body  contour  (which  Is  now  aligned  closely  w*‘ th 
the  horizontal  axis  In  the  zJB  k  space)  nearly  perpendicular  to  the 
image  of  the  centerline,  which  runs  between  h^  k  and  h2  jBk  (and  is 
still  a  s'  aight  line).  Accordingly,  a  simple  square  root  conformal 
transformaticn  may  be  applied  in  the  form 


zJC,k 


(z 


JB,k 


'  h2,JB,k) 


1/2 


(3.19) 


leading  to  the  hinge  point  alignment  shown  in  Figure  3.3.  Also  shown  in 
this  figure  is  the  resulting  body  surface  contour  for  the  simple  case  of 
a  sphere,  using  the  hinge  points  shown  in  Figure  3.2. 


Figure  3.3,  s  Plane  Hinge  Point  Images  and 
Body  Contour  (Sphere) 


Because  the  sequence  of  transformations  defi ned' above  is 


carried  out  independently  in  each  $  *  4>k  plane,  there  is  no  necessary 
correspondence  between  hinge  point  image  locations  in  these  planes, 
except  that  hg  jq  ^  *  0  and  Re  jc  k)  "  ®  ^or  ^  In  order  to 
minimize  the  discrepancies  that  must  arise  between  these  mappings  along 
the  centerline  (which  is  common  to  all  <j>^  planes),  a  final  stretching  is 
applied  in  each  plane  to  ensure  that  the  hinge  point  images  h-j  k 
coincide  in  the  =  £  +  in  space.  This  goal  is  attained  by  setting 

=  ak  2JC,k  *  (3.20) 

with  the  real  coefficients  a^  defined  by 

ak  a  hl,JC,l/hl,JC,k  *  (3‘21) 

This  simple  scaling  is  itself  a  conformal  transformation  and  thus 
preserves  the  orthogonal  nature  of  the  (£,n)  grid.  (Note  that  the  (£,n»6) 
space  is  not,  however,  necessarily  orthogonal.) 

It  is  important  to  note  that  while  the  final  images  of 
h.j  i  ^  and  hg  -|  ^  have  the  same  values  in  the  ^  space  that  only  those 
two  points  along  the  centerline  have  direct  correspondence  in  the  z-|  ^ 
space.  Because  of  different  scale  factors  that  arise  from  the  independent 
conformal  transformations  in  each  plane,  points  with  the  same  e 
value  do  not  necessarily  correspond  to  the  same  point  in  the  z-|  k  planes. 
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Although  the  conformal  mappings  in  each  plane  are  defined 
independently,  the  global  transformation  may  be  considered  continuous 
by  requiring  that  che  governing  parameters  of  the  transformation  be 
continuous  functions  of  4>  and  that  each  ^  plane  have  the  same  number 
of  hinge  points  (OC).  In  particular,  this  requires  that  the  functions 
a(<J>),  h2  j5 (4>) ,  hJ<+1  j(4>),  and  6^ (<p)  be  continuous. 

The  success  of  this  mapping  sequence  is  illustrated  in 
Figures  3.4  -  3.7.  Shown  in  these  figures  are  longitudinal  nosetip 
profiles  that  are  characteristics  of  low  altitude,  turbulent  ablation  of 
initially  spherical  carbonaceous  nosetips.  In  each  case,  the  hinge 
points  used  in  the  transformation  are  indicated  in  the  z-|  plane,  as  well 
as  the  body  contours  that  result  in  the  c  plane.  These  figures  indicate 
the  flexibility  inherent  in  this  conformal  mapping  procedure,  allowing 
any  arbitrary  nosetip  contour  to  be  mapped  onto  a  nearly  horizontal  line. 
Figures  3.5  and  3.6  represent  postulated  axisymmetric  nosetip  shapes  that 
have  been  tested  in  wind  tunnels:  the  Very  Mildly  Indented  Body  (VMIB), 
as  reported  by  Reeves,  Todisco,  Lin,  and  Pallone20  (Figure  3.5)  and  the 
PANT  Tri conic,  as  reported  by  Jackson  and  Baker21  f Figure  3.6).  Figure 
3.7  represents  a  profile  of  the  indented  nosetip  shown  in  the  Schlieren 
photograph  in  Figure  1.2,  which  was  recovered  from  a  flight  test. 

The  process  of  mapping  the  body  contour  onto  a  nearly  hori¬ 
zontal  line  is  relatively  insensitive  to  the  selection  of  hinge  points, 
as  long  as  the  hinge  points  approximate  the  body  shape  in  some  reasonable 
fashion.  Thus,  the  selection  of  hinge  points  is  easily  automated  by 
spacing  them  at  a  fixed  distance  along  inward  body  normals  (in  the  (x,y) 
plane)  from  body  points  equally  spaced  in  wetted  length. 
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3.3  TRANSFORMED  EQUATIONS  OF  MOTION 

Using  the  mapping  from  (x,y,4>,t)  space  to  (5»n*9,T)  space 
described  in  the  preceding  section,  the  governing  inviscid  equations 
(Equations  (3.1 )-(3.5)) may  be  transformed  to  (5.n>9,T)  coordinates  by 
application  of  the  chain  rule.  Recalling  the  functional  dependence  of 
the  transformation  defined  in  Equations  (3.11 )~(3.14),  the  appropriate 
chain  rules  take  the  forms 


3_  -  1_ 

3t  ‘  3t 


(3.22) 


3 _  k  3_  .  _  i_ 

3x  ■  ’X  35  nx  3q 


(3.23) 


3  -  -  1_  +  n  3 

3y  V35  ny  3n 


(3.24) 


a _  _ 
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r  <L_  +  n  3  .3 
S  35  n4>  3n  39 


(3.25) 


It.  is  convenient  to  define,  using  the  notation  of  Moretti17, 


g  «  «  Geiw  =  G(C  -  iS)  (3.26) 

3Z1 

and 

«  =  i0|Ol  =  4>1  +  i<|>2  (3.27) 

with 

(3.28) 


G  =  I  g  | 

«  =  arg(g) 


(3.29) 


C  -  cos  ( -to) 
S  a  sin  (-a)) 


From  the  definition  of  the  conformal  transformation,  it  follows  that 


g  =  ak 

JA 

^  9j  /2zJC,k 

(3.32) 

wi  th 

*  SJk 

(zj+l,k  "  1)/(zj,k  "  hj+l,j,k) 

(3.33) 

and 

$  * 

.  ’  tl 

akZJC,k  9 

JA 

(3.34) 

With  these  definitions,  the  partial  derivatives  required 

by  the  chain  rules  are  found  to  be 

£x  =  GC 

(3.35) 

«y  “  G~s 

(3.36) 

nx  -  -G's 

(3.37) 

n  =  GC 


(3.38) 


i'" 


Note  that  these  forms  verify  that  the  mapping  c  *  z  (z-j) 
in  any  $  =  constant  plane  is  conformal,  since  the  Cauchy-Riemann 
conditions 


5X  =  ny  (3.39) 

Cy  =  (3.40) 

are  satisfied. 

Circumferential  variations  of  the  mappings  are  accounted 
for  with  the  equations  (derived  from  a  Taylor  series  expansion) 


Z<p  =  £<j>  +  in<j)  *  "g(z2’“zl  /(4>2“4>1 ) 

9<j>  a  (92“gTg2  4  (Zg-z-j ))/(4>2“4>1 ) 

where 

C1  *  C(x1,y1,$1) 

Cg  =  £ (^2  »y2 j ^2 ) 

9-j  =  g(x1,y1,«1) 

§2 =  9(x2*^2,42^ 

with  ( x-j  ,yi  ,4>-j )  and  (xg »y2 »^2^  rePresentin9  computational  grid  points  in 
surrounding  4>  planes;  i.e.,  4>-j  =  <t>  -  A4> ,  d>2  =  4>  +  A<J>. 


(3.41) 

(3.42) 
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It  is  convenient  to  write  the  governing  equations  in  terms 
of  velocity  components  in  the  (C,n»9)  space.  Defining 

$  -  u?  +  vj  +  w£  (3.43) 

with  ?,  j,  and  k  being  unit  vectors  in  the  5,n,  and  0  directions, 
respectively,  the  new  velocity  components  may  be  written  in  terms  of 
the  cylindrical  velocity  components  (U,V,W)  as 

u  =  UC  +  VS  (3.44) 

v  =  -US  +  VC  (3.45) 

w  =  W  .  (3.46) 

In  terms  of  these  velocity  components,  the  governing  equations 
in  non-conservation  form  may  be  transformed,  using  the  chain  rules  de¬ 
fined  above,  to 

DP 

+  Y[G(u?  +  vn  +  v$2  -  u<t>-, ) 

+  (Sq  w£  +  j)  wn  +  w0  +  uS  +  vC)/y]  =  0  (3.47) 

^  +  vG(v$1  +  u$2)  +  vw(^4>2  +  n^-i  +  w0)/y 

-  Sw2/y  +  GpP^/p  =0  (3.48) 

{57  -  uG ( vd»1  +  u <t>2)  -  uw(^4>2  +  n^-,  +  u>0)/y 

-  Cw2/y  +  GpP^/p  =  0  (3.49) 
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(3.50) 


{*?  +  w(uS  +  vC)/y  +  p(^P5  +  r^Pn  +  P0)/py  =  0 


where 


D__  =  3 


UT 


f?  +  <Gu  + 


wVy) 

w/y  gg- 


3_ 

35 


+  (Gv  + 


wy*>  ?? 


i he  term  o)0  can  be  evaluated  from 
U)0  =  Im{g0/g) 
with  gQ  expressed  as 

9es9r9tS  • 


(3.51) 


(3.52) 


(3.53) 


3 COMPUTATIONAL  TRANSFORMATION 

Prior  to  obtaining  numerical  solutions  of  the  transformed 
governing  equations,  it  is  convenient  to  perform  an  additional  coordinate 
transformation  to  map  (5»n.0>T)  space  onto  a  rectangular  computational 
space  (X,Y,Z,T),  in  which  an  equally  spaced  mesh  can  easily  be  established 
to  facilitate  numerical  approximations  of  derivatives.  In  this  computa- 
•.  tonal  s*-  ■.  -tut,  the  coordinate  Z  is  selected  so  as  to  be  0  on  the  body 
surface  and  1  on  the  outer  boundary  (bow  shock  wave)  of  the  region  of 
interest.  S-’.'rilarly,  Y  is  defined  as  being  0  on  the  centerline  and  1  at 
the  downstr  boundary  of  the  region  to  be  computed.  X  is  directly 
proportional  to  the  circumferential  coordinate  8. 
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This  computational  transformation  is  described  mathematically 


as 

X  =  8/2ir  (3.54) 

Y  =  S/5l(9)  (3.55) 

z  =  [n-b(5,e)]/[c(5,e.T)-b(e,e)]  (3.56) 

T  =  t  (3.57) 

where  the  body  surface  is  described  as 

n  «  b(e,e)  (3.58) 

and  the  bow  shock  position  as 

n  «  c(C,0,t)  .  (3.59) 


The  downstream  boundary  is  defined  by 


5  =  eL(e) 


(3.60) 


Because  the  position  of  the  bow  shock  varies  with  time 
during  the  solution  of  the  time-dependent  problem,  the  computational 
grid  also  varies  with  time,  but  always  maintains  equally  spaced  points 
between  the  body  and  bow  shock  (in  Z). 

To  transform  the  governing  equations  into  the  computational 
coordinates,  the  following  chain  rules  are  applied: 


9_ 

3t 


(3.61) 


9__  Y  3 

35  5  W 


(3.62) 
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where 


9_ 

an 


3__ 

90 


from 


s  Z  £_ 

4n  3 1 

(3.63) 

=  y  l_+v  —  +  7  1_ 

A0  3X  T0  3Y  0  3Z 

(3.64) 

*  1/2tt 

a  i/SL(e) 

*  'VYrC, 

5  l9 

•  V[c(s,0fT)-b(c,e)] 

*  -2nt(l-Z)b,  +  Zc?] 

-  -^[(l-Z)b9  ♦  ZCq) 

=  -ZZnS 

and  shock  slopes  in  the  transformed  s'ace 

being  determined 

*  (cybx  -  s)/(sybx  +  c) 

(3.65) 

*  sy^f^bx +  -  Vc +  "♦ 

(3.66) 

■  <^sx  -  SV(Sysx  +  C) 

(3.67) 

■  *,*/&„ +  e>  ■  Vs +  r>$ 

(3.68) 

Using  Equations  (3.54)-(3.57)  to  define  the  computational 
transformation  and  selecting  equal  mesh  intervals  in  X,Y,  and  Z  leads 
to  points  equally  spaced  in  r\  (between  body  and  shock)*  in  £  (along 
the  body),  and  in  u  (circumferentially).  It  is  possible,  however,  to 
use  stretching  functions  in  the  definitions  of  X,Y,  and  Z  in  order  to 
concentrate  mesh  points  in  certain  regions  if  so  desired,  while  main¬ 
taining  equal  spacing  in  X,Y,  and  Z.  For  example,  Moretti17  and 
Kutler,  Chakravarthy,  and  Lombard18  have  used  such  stretching  functions 
to  concentrate  grid  points  near  the  body  surface  to  facilitate  viscous 
calculations.  For  the  nosetip  inviscid  flow  problem,  such  stretching 
is  not  deemed  necessary,  and  the  simpler  definitions  of  X,Y,  and  Z,  as 
presented  above,  are  used. 

The  governing  equations  written  in  the  computational  space 
may  now  be  expressed  as 

{jf  +  YG[Y^uy  +  Z?uz  +  Znvz  +  E  +  (uS  +  vC)/Gy] 

+  T^<j>(Y5wY  +  ZFWZ^  +  n4>ZnwZ  +  X0WX  +  Y0wY 
+  ZgW-z3/y  =  0  (3.69) 

+  vGD  +  vwF/y  -  Sw2/y  +  GpfY^Py  +  l^z)/P  =  0  (3.70) 

-  uGD  -  uwF/y  -  Cw2/y  +  GpZ nPz/p  =  0  (3.71) 

^+w(>,S  +  vC)/y  +  p[(Z^  +  Znnt  +  Z9)Pz 

+  (^  +  Y3)Py  +  XePx)/py  »  0  (3.72) 
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with 


Ds  -  n 

ET  “  0 


D _  _  3 _  .  a  3  +  o^  +  r  9 

DT'FT+A3l+B37+C9l 

A  =  ZT  +  (Gu  +  wC^/yjZ^  +  (Gv  +  wn(J>/y)Zn  +  wZQ/y 
B  ■  (Gu  +  w^/yjY^  +  wYQ/y 


C  =  wXQ/y 


(3.|73) 


D  *  v<J>,  +  u<}>2 
B  =  v<j>2  "  u*1 

F  =  V2  +  Vl  +  “9  * 

Typical  grids  (in  physical  space,  <J>  =  constant)  that  result 
from  this  computational  transformation  are  shown  in  Figures  3.8-3.10. 

(In  these  figures,  the  bow  shock  shape  used  in  defining  the  computational 
region  is  an  assumed  initial  bow  shock  shape.)  For  clarity  in  these 
figures,  a  coarser  grid  is  shown  than  would  actually  be  used  in  the 
calculation  of  a  flow  field  about  such  bodies. 

It  is  important  to  note  in  these  figures  that  the  £  =  constant 
lines  are  indeed  nearly  normal  to  the  body  surface,  as  is  expected  when 
the  image  of  the  body  contour  is  nearly  horizontal  (n  =  constant)  and 
the  mapping  is  conformal.  The  generation  of  such  grids  was  the  goal 
in  the  development  of  the  mapping  function  presented  here,  and  will 
greatly  expand  the  range  of  nosetip  shapes  that  can  be  successfully  computed. 
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3.5 


CHARACTERISTIC  RELATIONS 


The  numerical  procedures  to  be  used  at  body  and  bow  shock 
points  (which  will  be  discussed  in  Sections  5.3  and  5.4)  are  based 
on  the  characteristic  compatibility  relations  resulting  from  the 
governing  system  of  partial  differential  equations,  Equations  (3.69)- 
(3.73).  Accordingly,  the  appropriate  forms  of  these  compatibility 
relations  at  the  body  and  bow  shock  are  derived  in  this  section. 

In  the  theory  of  partial  differential  equations  a  charac¬ 
teristic  surface  is  a  surface  across  which  the  derivatives  of  the 
dependent  variables  may  be  indeterminant  in  the  direction  normal  to  the 
surface.  The  characteristic  compatibility  condition  is  a  linear  combi¬ 
nation  of  the  governing  equations  valid  along  this  characteristic  surface. 

For  use  in  this  analysis,  characteristics  in  the  (Z,T) 
reference  plane  are  of  interest;  X  and  Y  derivatives  appearing  in  the 
governing  equations  will  be  treated  as  forcing  functions.  Reduction  of 
the  four-dimensional  (X,Y,Z,T)  problem  to  two  dimensions  results  in  a 
characteristic  curve,  rather  than  a  characteristic  surface. 

The  governing  equations  may  be  rewritten  as 

PT  +  APZ  +  tG(Z£uZ  +  ZnvZ^ 


+  +  n(j)Zn  +  ZQ)wz/y  =  R1 

(3.74) 

+  Auz  +  GpZ^Pz/p  =  R2 

(3.75) 

+  Avz  +  GpZnPz/p  =  R3 

(3.76) 

+  Awz  +  p(^Z?  +  rt(j)ZT1  +  Z0)Pz/py  =  R4 

(3.77) 
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where 


R1  =  “ [BP Y  +  CPX  +  yG(Y^uy  +  E  +  V/Gy) 

+  +  XeWX  ••  YgWy)/y] 

Rg  =  -[Buy  +  Cu^  +  vGD  +  vwF/y  -  Sw2/y 
+  GpY^/p] 

Rj  =  -[Bvy  +  Cvx  -  uGD  -  uwF/y  -  Cw2/y] 

R4  =  -[BwY  +  Cw^  +  w(uS  +  vC)/y 
+  p{(^Y?  +  Yq)Py  +  X0Px}/py] 

(The  equation  for  entropy  convection.  Equation  (3.73),  is  not 
considered  here,  since  it  is  known  that  the  characteristic  resulting 
from  its  inclusion  is  simply  a  streamline.  While  a  streamline  is  a 
valid  characteristic,  it  is  not  of  immediate  interest  for  this  application.) 

Defining  the  characteristic  curve  as 

f(T,Z)  =  0  (3.78) 

the  normal  to  this  curve  is 

ft  «  Vf  =  (fT,fz)  (3.79) 

and  the  characteristic  slope  may  then  be  defined  as 

X  =  -  P-  .  (3.80) 

rZ 
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The  characteristic  compatibility  condition  is  written  as  a  linear 
combination  of  the  governing  equations,  where  JLj ,  Jt^,  *3  and  *4  are 
the  as  yet  undetermined  multipliers  for  Equations  ( 3 . 74) - ( 3 . 77 ) , 
respectively.  Combining  terms,  the  compatibility  condition  can  then 
be  written  as 

Vt  +  W1A  +  gP(12zc  +  +  Wc  +  Vn  +  Z8)/w1,,2 

+  ^2ux  +  (^ySZ^  +  S-gA^z  +  AgVj  +  (^ySZ^  +  ^3A)vz 

4 

+  Vt  +  wiy(V5  +  Vn  +  ze)/y  +  *4A]wz  =  ?  *iRi  *  (3-81) 


The  terms  involving  derivatives  of  P  may  be  regarded  as  a 
directional  derivative  in  the  direction  ftp  where 

*  [^i»{^iA  +  GpCi^Z^  +  ^Z^/p  +  fc4p(  Vs 

+  Vn  +  V/py}]  •  (3.82) 

Similarly,  derivatives  of  u,v,  and  w  may  also  be  viewed  as  directional 
derivatives  in  the  directions  ancl  V 

tf2  =  I^.^YGZ^  +  A2A1  (3.83) 

i?3  =  [  WGZn  +  *3A]  (3.84) 

^4  =  t V^l  y(^4>Z?  +  r'4)Zn  +  Z9^y  +  ^A)]  *  (3.85) 
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For  Equation  (3.81)  to  be  valid  along  the  characteristic, 
these  directional  derivatives  must  not  have  any  component  along  the 
direction  of  the  normal  to  the  characteristic  curve  (in  which  direction 
the  derivatives  may  be  indeterminant).  These  conditions  may  be  expressed 
as 


ft  •  ftj  -  ft  •  &2  -  ft  •  &3  *  ft  •  ft4  «  0 


(3.86) 


g  g 

Noting  that  gY  a  37  »  this  system  of  equations  may  be  written  in 
matrix  form  as 


A-X 


YGZ. 


yG7. 


GpZ^/p  GpZn/p  P(^Ze+T^Zn+Z0)/py 


A-X 


Y(V€+Vn+V/y  0 


A-X 


A-X 


\ 


ri 


'3 

^4  / 


=  ft 


(3.87) 


For  a  solution  to  this  system  of  homogeneous  equations  to 
exist,  it  is  necessary  that  the  determinant  of  the  coefficient  matrix 
vanish.  Furthermore,  any  one  of  the  four  unknowns  may  be  scaled  arbi¬ 
trarily.  Expansion  of  the  determinant  results  in  the  following  algebraic 
equation: 


(A-X)4  [(A-X)2  -  a2G2Z52  -  &2G2Zn2 

-a!(5+Z;  +  Vn  +  V2/*2'*  0 


(3.88) 


■ - r’rTFVB-? 


where  the  isentropic  speed  of  sound  is  defined  from 

a  =  (yp/p)1/2  .  (3.89) 

The  four  roots  to  this  equation  are 

X  =  A  (redundant  root)  (3.90) 

and 

*  ■  A*  a  I5!(V  +  Zn!>  *  p  (5„ZE  +  *  Ze)2]1/2  .  (3.91) 

The  redundant  root  X  =  A  simply  shows  that  streamlines  are  characteristic 
directions,  but,  as  stated  earlier,  this  relation  is  not  of  immediate 
interest.  Thus,  the  characteristic  slopes  being  sought  are  those  defined 
by  Equation  (3.91). 

To  evaluate  the  unknown  multipliers  it  is  convenient  to 
select  =1;  it  then  follows  that 

l2  «  yGZ?/(X-A)  (3.92) 

5-3  =  yGZn/(X-A)  (3.93) 

+  r^Zn  +  Z0)/(X-A)y  .  (3.94) 

The  compatibility  condition  will  then  take  the  final  form 
Py  +  XPZ  +  %2  (uy  +  *uz)  +  53(v-j.  +  Xvz) 

4 

+  Mwt  +  Xw7)  =  E  .  (3.95) 

4  1  L  i=l  1  1 
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To  derive  forms  of  this  relation  valid  at  the  body,  it  is 
first  necessary  to  write  the  kinematic  boundary  condition  in  (£>n»9) 
space.  Denoting  the  body  normal  as 

nb  =  -Gb^i  +  Gj  +  (T^-S^-bgJ/y  k  (3.96) 

the  boundary  condition  becomes 

-Gub^  +  Gv  +  Wn^-C^-b^/y  =  0  .  (3.97) 

The  coefficient  A,  defined  as 

A  =  ZT  +  (Gu  +  w^/y)Z?  +  (Gv  +  wn({)/y)Zn  +  wZQ/y 

can  be  shown  to  vanish  at  the  body  since,  with  Z  =  0, 


Z5  ■  "Znb5 


Ze  ■  'znbe 


and  thus 

A  *  Z^-Gub^  +  Gv  +  w(n(f)-?(j>b?-be)/y]  £  0 

from  Equation  (3.97).  Choosing  X  <  0  at  the  body  and  simplifying  the 
expression  for  X  yields 

*b  -  -aZn[G2(Hb52)  +  (TV9=E-b9)W!1/2 


(3.98) 


and  the  compatibility  condition  becomes 

PT  +  XPZ  +  YZnt-Gb^uz  +  Gvz  +  (n^-^bg 

4 

-b0)wz/y]  =  E  (3.99) 

since  J^Uj  +  ^3Vj  +  t^Wy  =  0  from  the  time  invariant  boundary  condition, 
Equation  (3.97). 

At  the  shock,  with  X  >  0  and  Z  =  1 ,  it  follows  that 

xs  *  A  +  iln  tG2(l+c5!)  +  (TyS^-Cg)2//]1'2  (3.100) 

since 

h  =  'ZncE 

z9  *  ‘Ve 

and  A  no  longer  vanishes.  The  appropriate  compatibility  condition  at  the 
shock  is  given  in  the  general  form  of  Equation  (3.95). 

Specific  application  of  these  characteristic  relations  for 
boundary  point  calculations  will  be  presented  in  Section  5.0. 

3.6  TREATMENT  OF  THE  SINGULAR  CENTERLINE 

In  the  (x,y,4>)  cylindrical  coordinate  system  the  x-axis  (y  -  0) 
is  a  singular  line,  where  the  coordinate  ef  is  multi-valued.  Along  this 
axis  the  governing  Equations  (3.69)-(3.73)  must  take  on  different  forms 
valid  along  this  axis,  eliminating  indeterminate  terms  that  result  from 
the  singularity.  The  modified  governing  equations  that  then  result 
involve  second  derivatives  of  the  dependent  variables,  such  as  P^y»  etc. 
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In  order  to  avoid  these  second  derivatives,  other  approaches 
to  the  centerline  problem  may  be  used.  For  example,  the  governing  equa¬ 
tions  can  be  reformulated  in  a  local  Cartesian  coordinate  system,  which 
does  not  exhibit  singular  behavior.  Approximation  of  derivatives  in 
the  Cartesian  space  would,  however,  require  extensive  interpolation  on 
the  data  at  the  computational  grid  points  which  are  not  aligned  with  the 
Cartesian  coordinates. 

In  this  analysis,  a  set  of  governing  equations  based  on  the 
Cartesian  approach  at  the  centerline  is  developed  which  minimizes  the 
need  for  interpolation,  while  simultaneously  avoiding  the  approximation 
of  second  derivatives.  (For  the  three-dimensional  conformal  mapping 
approach  developed  in  this  research  effort,  the  approximation  of  second 
derivatives  at  the  centerline  is  made  particularly  difficult  by  the  fact 
that  the  transformations  used  in  each  $  plane  are  independent  and  thus 
are  not  continuous  across  the  centerline.) 

To  develop  the  form  of  the  governing  equations  desired  at 
the  centerline,  consider  a  Cartesian  reference  frame  (x^,x2,x3),  oriented 
with  the  (x,y,4>)  cylindrical  system  as  shown  in  Figure  3.11,  and  defined 
by 


=  X 

(3.101) 

=  y  cos  <b 

(3.102) 

=  y  sin$ 

(3.103) 
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Derivatives  in  the  Cartesian  frame  may  now  be  expressed  in  terms  of  the 

cylindrical  frame,  without  any  indefinite  forms  appearing,  by  carefully 

selecting  the  values  of  <j>  for  which  certain  derivatives  are  evaluated. 

For  example,  all  derivatives  |y~  are  evaluated  as  cos$  ^y  in  the  4>  «  0 

^3  3 

and  i>  -  tt  planes,  and  all  derivatives  |y-  are  evaluated  as  sin<j>  in 

3 

the  4>  =  it/2  and  <j>  -  3tt/2  planes.  These  simple  forms  result  since 


lim  1  3  _ 


y-K)  y  34>  .  3y3  j) 


(3.107) 


which  has  a  finite  value  (if  and  are  bounded,  as  is  implicitly 


3x. 


assumed  in  this  analysis). 

Starting  with  the  governing  inviscid  equations  written  in  a 
Cartesian  coordinate  system  and  applying  these  chain  rules,  a  system  of 
equations  in  cylindrical  coordinates  results  that  is  valid  along  the 
centerline  and  does  not  involve  any  second  derivatives.  The  resulting 
equations  do,  however,  have  some  terms  that  must  be  evaluated  in  the 
<|>  =  0  or  4>  *  it  planes,  and  others  that  must  be  evaluated  in  the  4>  =  tt/2 
or  <j>  =  3ir/2  planes.  (Because  of  this  form  of  the  equations,  it  is 
necessary  in  the  numerical  solution  to  require  a  computational  grid  that 
includes  these  four  4)  planes.) 

Transforming  these  special  onnations  in  cylindrical  coordi¬ 
nates  tc  (C,ti,0,t)  space  and  then  to  the  (X.Y.Z.T)  computational  space, 
and  writing  the  equations  in  terms  of  the  transformed  velocity  components 
results  in  the  final  forms  of 


[P,,  +  A^z  +  B1Py  +  yG  Ucuz  +  Y^Uy  +  Znv7  +  *  0,* 

-v  [Gu(ZpPz  +  Y  PY)  +  yG(Z5uz  +  Y^Uy  +  V4>2)]4*  |L  s  0 


(3.108) 
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(3.109) 


[{Uy  +  A-jU^  +  B^Uy  *  vGD  +  +  Y^Py)/p)cos^]^  B 

+  t-Gu(Z^wz  +  Y^Wy)sin(j)](j)  _  u  3tt  »  0 


[vT  +  A-jVz  +  BlVy  -  uGD  +  pGZnPz/p](()  .  Q ^ 

+  [Gu(Z-v7  +  YPvv  -  u<J>9)].  _  tt  3tt  =  0 
5  i  5  Y  2 '  <t)  - 

[{w-j-  +  A^wz  +  B-jWyJcosif)]^  .  Q|  +  [(Gu(Z^uz  +  Y^Uy  +  v^) 
+  pG(Z^Pz  +  Y^Py )/p)si r»4>] ^  _  it  3tt  =  0 


(3.110) 


(3.111) 


CsT  +  A1SZ  +  BlSY1<j)  *  0,ir  *  [Gu(Vz  +  YCSY^4>  =  |  ,  Q  =  0  (3.112) 


where 

A1  =  Zt  +  G(uZ^  +  vZn) 

B1  =  GuY? 

The  characteristic  compatibility  conditions  required  at  the  body  and 
shock  points  on  the  centerline  may  be  formed  as  a  linear  combination  of 
these  special  centerline  equations  following  the  same  procedure  presented 
in  Section  3.5  for  regular  points.  Special  forms  of  the  characteristic 
slopes  X  and  the  multipliers  may  be  derived  at  the  centerline  as 
follows. 
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Consider  the  unit  body,  normal  at  the  centerline.  In  general 


the  unit  body  normal  at  points  off  the  centerline  may  be  written  as 
a  a  -Sbgi  +  G3  +  (y6»bg-b9)/y  fc 

n  *  [6* (1+b^2 )  +  (ry^-b^Vy2)1'2  '  (3.113) 

At  the  centerline,  of  course,  a  different  expression  must  be  used  for 

A 

the  k  component  of  this  vector.  But  now  consider  this  unit  vector  at 

AAA 

the  centerline  in  the  $  a  0  plane  (with  unit  vectors  i.,j.j,k,)  and  in 

AAA 

the  <|>  *  tt/2  plane  (with  unit  vectors  i2,j2,k2),  as  depicted  in  Figure 
3.12.  This  unit  body  normal  can  then  be  written  in  the  two  equivalent 
forms 


(3.114) 

(3.115) 

where  D-j  and  D2  are  the  respective  denominators.  But  noting  that 

and  j2  are  coincident  with  the  centerline  and  that  Equations  (3.114) 
and  (3.115)  represent  the  same  vector,  it  follows  that 


-G-jb^i-j  +  G-jj-j  +  (n^-E^b^  -  b@^)/y^  k^ 

- - d; 

~  -G2b£2i2  +  G2j2  +  (n^i-E^ibg^  -  b0^)/y^  k.| 


(3.116) 
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where  A  at  $  =0  is  now  given  by 


A  *  A1  -  GwZnC£2  .  (3.123) 

The  multipliers  at  the  body  become 

£2  *  -yGZ^/X  (3.124) 

=  yGZ  /X  (3.125) 

<5  n 

t4  =  -yGZ^/X  (3.126) 

and  at  the  shock 

Jt2  =  -yGZtic^1/(X-A)  (3.127) 

*3  =  YGZn/(X-A)  (3.128) 

t4  =  -YGZnc?2/(x-A)  .  (3.129) 


While  the  expressions  presented  above  have  been  derived 
assuming  that  (  )^  refers  to  <J>  *  0  and  (  )2  refers  to  <J>  ■-  tt/2,  these 

forms  are  equally  valid  for  (  )^  representing  <j>  =  n  and  (  )2 
representing  4>  =  3tr/2. 
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SECTION  4 

CALCULATION  OF  EMBEDDED  SHOCKS 

The  method  selected  in  this  research  for  the  calculation 
of  embedded  shocks  is  the  shock-capturing  approach,  in  which  shock 
waves  (and  other  discontinuities,  such  as  slip  lines)  are  computed 
automatically,  albeit  approximately.  Two  methods  of  shock-capturing 
are  examined  in  this  section:  the  conservation  (shock-smearing) 
approach  and  the  X-differencing  (non-conservation)  approach.  The 
relative  merits  of  these  two  methods  are  compared  by  developing 
axisymmetric  versions  of  both  procedures,  described  in  Sections  4.1 
and  4.2,  and  assessing  the  abilities  of  each  scheme  to  compute  in- 
viscid  shock  layer  flows  with  embedded  shocks.  Section  4.3  details 
the  comparisons  of  these  calculations  to  experimental  data,  which 
show  the  X-differencing  approach  to  be  superior.  Finally,  in  Section 
4.4,  the  X-differencing  scheme  is  extended  to  three  dimensions. 

4.1  CONSERVATION  LAW  APPROACH  TO  SHOCK-CAPTURING 

The  theory  behind  the  conservation  law  approach  is  to 
reformulate  the  governing  partial  differential  equations  in  terms 
of  dependent  variables  that  appear  naturally  in  the  integral  con¬ 
servation  laws.  The  resulting  dependent  variables  then  represent 
quantities  that  are  reminiscent  of  the  quantities  conserved  across  a 
discontinuity  from  the  Rankine-Hugoniot  conditions  (mass,  momentum, 
and  energy  flux).  Hopefully,  these  new  dependent  variables  will  be 
continuous  across  the  discontinuity,  and  thus  a  numerical  solution 
can  be  obtained  directly  without  special  treatment  for  the  discontinuity. 
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It  must  be  noted,  however,  that  the  dependent  variables  in 
the  differential  conservation  formulation  are  strictly  continuous  only 
if  the  discontinuity  is  perfectly  aligned  with  the  coordinate  mesh 
used  in  the  calculations  and  if  the  discontinuity  is  stationary.  There¬ 
fore,  for  shocks  that  are  not  aligned  with  the  mesh  or  that  are  moving 
(such  as  during  the  transient  phase  of  a  time-a-ymptotic  calculation), 
the  dependent  variables  are  not  continuous  and  the  conservation  form 
of  the  governing  differential  equations  is  not  strictly  valid. 

To  illustrate  these  points,  consider  a  stationary  shock 
inclined  at  an  angle  o  to  a  two-dimensional  Cartesian  mesh,  as  shown 
in  Figure  4.1.  The  conservation  form  of  the  steady  inviscid  continuity 
equation  is 

(pu)x  +  (pv)y  »  0  (4.1) 

where  u  and  v  are  the  x-  and  y-  velocity  components,  respectively. 


Figure  4.1.  Steady  Inclined  Shock 
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(4.2) 


The  Rankine-Hugoniot  conditions  for  this  case  require  that 
P^l  *  p2°2 

v1  =  Vg  (4.3) 

where  (  )-j  denotes  the  low  pressure  side  of  the  shock,  and  (  )2 

denotes  the  high  pressure  side.  Since 

u  -  u  cos  a  -  v  sin  o  (4.4) 

v  s  u  sin  o  +  v  cos  o  (4.5) 

the  continuity  relation  given  by  Equation  (4.2)  becomes 

p-| u.j  =  p2u2  +  tan  o  (p-jV-j  -  p2v2)  .  (4.6) 

Clearly,  the  conservation  variable  pu  will  be  continuous  across  the 
shock  only  if  a  =  0  (i.e.,  if  the  shock  is  aligned  with  the  coordinate 
system). 

Similarly,  consider  the  case  of  a  normal  shock  (o  =  u/2) 
moving  to  the  left  with  velocity  W,  as  illustrated  in  Figure  4.2. 

Equation  (4.2)  becomes 

p.j(u.|  +  W)  *  p2(u,>  +  ’  (4-?) 

Hence  the  conservation  variable  pu  will  be  continuous  across  the  shock 
only  if  the  shock  velocity  were  to  vanish. 

A  discussion  of  the  dependence  of  conservation  shock¬ 
capturing  results  on  the  orientation  of  the  discontinuity  relative  to 
the  coordinate  mesh  may  be  found  in  MacCormack  and  Paullay25. 


w 


Figure  4.2.  Unsteady  Normal  Shock 

Despite  the  lack  of  continuity  of  the  conservation  variables 
across  a  shock  in  the  general  case,  however,  it  can  reasonably  be 
expected  that  the  conservation  variables  will  be  "smoother"  than 
the  primitive  variables  (p,p,etc.)  across  a  shock.  Thus,  conservation 
form  calculations  may  have  the  potential  of  "automatically"  computing 
shocks  in  cases  where  the  non-conservation  (Eulerian)  formulation, 
wi+hout  * -differencing  or  shock-fitting,  would  fail. 

The  calculation  of  discontinuities  with  the  conservation 
formulation  smears  the  discontinuities  over  several  mesh  intervals 
and  also  introduces  oscillations  into  the  calculation  at  discontinuities. 
The  conservation  approach  must  then  be  viewed  as  an  approximate  method 
for  computing  embedded  shocks  since  the  discrete  discontinuity  is  smeared 
out;  the  results  obtained  with  this  approach  will  thus  be  mesh  dependent. 
This  approach  requires  a  fine  computational  mesh  to  obtain  accurate 
approximations  to  embedded  shocks. 
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The  presence  of  spurious  oscillations  in  the  conservation 
calculations  may  require  special  numerical  treatment  to  avoid  failure 
of  the  computation.  The  procedures  used  in  this  analysis  to  control 
these  oscillations  are  detailed  later  in  this  section. 

The  conservation  form  of  the  governing  equations  can  be 
derived  in  many  forms;  the  recommended  formulation  for  the  calculation 
of  embedded  shocks  is  the  "strong"  conservation  form,  in  which  no 
undifferentiated  terms  appear,  leading  to  the  overall  conservation  of 
mass,  momentum,  and  energy,  as  discussed  by  Vinokur26.  (Note,  however, 

that  this  is  a  global  conservation  of  mass,  momentum,  and  energy.)  The 

» 

totally  "strong"  conservation  form  cannot  be  obtained  for  the  axi sym¬ 
metric  equations,  however. 

The  axisymmetric  conservation  equations  in  the  computational 
coordinate  system  may  be  written  as 

^T  +  +  +  ^  =  ^  (4.8) 

where  the  vector  quantities  are  defined  as 

’  p  ' 

PU 

< 

pv 

k  6  4 

•  pA  ' 

pUA  +  G(CZ-  -  SZ  )p 

(  _  _  n  ► 

pVA  +  G(SZ?  +  CZn)p 

>  eA  +  G(uZ^  +  vZ^Jp 


t  = 


*  _  1 
®  "  ¥T_ 
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pu  ' 
Cp  +  pull 
Sp  +  puV 
u(e  +  p) 


^  e  +  p  J 

with  e  representing  the  total  energy  per  unit  volume, 

e  =  pH  =  p [h  +  i  (U2  +  V2)]  (4.9) 

and  the  computational  coordinate  Y  being  redefined  for  the  axisymnetric 
case  as 

Y  =  £  (4.10) 

and  the  coefficient  A  representing 

A  =  Z  +  GuZr  +  GvZ  .  (4.11) 

x  5  n 

After  each  computational  step  (predictor  or  corrector),  the 
conservation  vector  F  must  be  decoded  to  recover  the  primitive  flow 
variables.  The  quantities  p,u»v,  and  h  can  be  determined  directly  f^om 
the  pressure  p  a. id  entropy  s  can  be  found  directly  from  these 


quantities  for  an  ideal  gas,  while  an  iterative  decoding  procedure  may 
be  required  for  equilibrium  real  gas  thermodynamics. 

The  use  of  the  conservation  formulation  of  the  governing 
equations  can  lead  to  the  presence  of  spurious  oscillations  in  the 
vicinity  of  a  smeared-out  discontinuity.  Typically,  the  magnitude  of 
the  oscillations  tends  to  increase  as  the  shock  strength  increases, 
and  if  the  oscillations  are  undamped,  the  calculation  can  quickly  fail 
as  the  oscillations  spread  throughout  the  shock  layer  being  computed. 

Most  conservation  law  techniques  use  some  form  of  numerical 
damping  (either  implicit  or  explicit)  to  control  these  oscillations 
(which  are  mesh  dependent).  In  this  analysis,  a  simple  damping 
technique  has  been  used,  which  is  shown  to  be  equivalent  to  an  arti¬ 
ficial  viscosity,  similar  to  that  used  in  other  conservation  formula¬ 
tions,  such  as  by  Lax  and  Wendroff 21 . 

To  illustrate  the  damping  procedure  used,  consider  the 
simple  hyperbolic  equation 


f .  +  u  f  +  v  f  *0 
t  ox  o  y 


(4.12) 


The  numerical  solution  to  this  equation  is  first  advanced  one  time 

step  using  the  MacCormack  7  predictor-corrector  finite  difference 

scheme,  described  in  Section  5.6.  Following  the  corrector  stage,  a 

k+1 

weighted  averaging  of  the  solution  f  is  performed,  as 


2k+l 


-k+1 


k+1 


rk+1 


r  1  =  ar  1  +  6[f  +  f  i  ] 

nm  nm  p  1  n+1  ,m  n-1  ,mJ 


+  Y^n+L.i  +  f„+i  il 

n,m+i  n,m- 1 


(4.13) 
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where  fnm  =  f(tQ  +  kAt,  xQ  +  nAx,  yQ  +  mAy)  and  subject  to  the  con¬ 
straint  that 

a  +  2$  +  2y  B  1  .  (4.14) 

Application  of  this  averaging  procedure  to  the  numerical 
solution  of  Equation  (4.12)  can  be  shown  through  Taylor  series  expansion 
to  be  equivalent  to  the  solution  of  the  modified  equation 

ft  +  uofx  +  vofy  *  8(ix)Vit  fxx  +  Y(ay)Vit  fyy 

+  0(At2,  Ax2,  Ay2)  .  (4.15) 

The  coefficients  of  the  second  order  terms  are  thus  similar  to  viscosities. 
For  consistency  of  these  viscous-like  coefficients,  the  coefficient  y 
can  be  selected  to  be 

y  =  B(Ax/Ay)2  .  (4.16) 

To  more  closely  simulate  physical  viscosity,  this  damping  is  applied 
only  to  the  two  momentum  equations  appearing  in  the  ax i symmetric  con¬ 
servation  system.  Equation  (4.8). 

While  the  damping  formulation  described  above  can  be  helpful 
in  controlling  the  oscillations  that  arise  in  the  calculation  of  embedded 
shocks  with  the  conservation  equations,  it  cannot  completely  compensate 
for  large  oscillations.  If  the  oscillations  are  severe  enough,  the 
calculation  will  typically  fail  when  a  negative  pressure  (p<0)  is 
encountered  in  the  decoding  of  the  conservation  variables.  Since  this 
condition  might  occur  only  during  the  transient  phase  of  a  time-dependent 
calculation,  and  not  in  the  steady-state,  a  second  damping  (or  smoothing) 
technique  is  applied  locally  when  required  to  eliminate  the  p<0 
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condition  and  allow  the  calculation  to  continue.  At  any  point  Y  *  mAY 
and  Z  *  nAZ  where  p<  0  is  encountered,  the  entire  F  vector  at  that  Y 
location  is  smoothed  in  Z,  from 

fnm'T  tTn+l,m*  2?n,m  + 

This  is  the  same  smoothing  function  used  by  Solomon,  et  al.12  when  the 
p<0  condition  is  encountered  in  the  conservation  calculation  of  embedded 
shocks  in  steady  supersonic  flow  on  reentry  vehicle  afterbodies.  This 
local  smoothing  function  reduces  the  numerical  scheme  to  first  order 
accuracy  in  AZ,  but  has  been  found  helpful  in  overcoming  transient 
difficulties  in  the  time -dependent  calculation. 

4.2  X-DIFFERENCING  APPROACH  TO  SHOCK-CAPTURING 

The  x-differencing  approach  to  shock-capturing  is  fundamentally 
different  from  the  conservation  approach  in  that  its  shock-capturing 
ability  comes  not  from  a  special  formulation  of  the  dependent  variables, 
but  rather  from  careful  treatment  of  the  approximations  to  spatial 
derivatives.  By  constructing  finite  difference  approximations  that 
accurately  model  the  domain  of  dependence  of  each  point  being  computed, 
the  resulting  finite  difference  equations  admit  "discontinuous"  solu¬ 
tions  that  approximate  the  structure  of  physics'!  solutions  with  embedded 
shocks. 

The  X-scheme  is  formulated  in  terms  of  the  non-conservation 
governing  equations,  in  which  the  dependent  variables,  in  two- di mens i ons , 
are  P,u,v,  and  s.  To  illustrate  construction  of  the  X-scheme,  consider 
a  one-dimensional,  time-dependent  inviscid  flow,  whose  governing  equations 


are 


Pt  *  uPx  +  YUx  *  0 


(4.17) 


ut  +  uux  +  p/pPx  *  0  (4.18) 

st  +  usx  *  0  *  (4.19) 

This  system  of  equations  has  two  characteristic  directions  (aside  from 
the  particle  path,  which  is  not  of  interest  in  this  application),  de¬ 
fined  as 

X-|  =  u  -  a  (4.20) 

X2  -  u+a  .  (4.21) 

Note  that  supersonic  flow  implies  X-j  >  0  and  X2>0  and  that  subsonic 
flow  produces  X-j  <  0  and  X2>0,  as  illustrated  in  Figure  4.3.  The  signs 
of  these  characteristic  slopes  at  any  point  deF’ne  the  direction  of 
the  domain  of  dependence  at  that  point. 

The  basis  of  the  X-scheme  is  to  rewrite  Equations  (4.17) 
and  (4.18)  in  such  a  manner  that  the  domain  of  dependence  information 
inherent  in  the  characteristic  slopes  can  easily  be  incorporated  into 
the  resulting  finite  difference  scheme.  To  this  end,  the  governing 
equations  may  be  written  as 

Pt  +  ?  +  ^2Px2^  *  y(X2Ux2  “  ^i^X] )/2a  =  0  (4.22) 

u^  +  t;-  (X^ ux-j  +  X2Ux2)  a(X2Px2  ~  ^lPxi ~  (4.23) 

which  are  entirely  equivalent  to  Equations  (4.17)  and  (4.18),  with 
Pxi  =  ^X2  anc*  uxi  =  UX2" 
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X 

a.)  SUBSONIC  FLOW 


x 

b.)  SUPERSONIC  FLOW 

Figure  *1.3.  Character  uric  * 

Slopes  in  x-t /  Plane 

I 

In  a  finite  difference  representation  of  these  equations, 
PX1  and  Pxg  (or  and  ux2)  need  not  be  the  same  approximation  to  Px 
(or  ux).  By  constructing  second-order  accurate  one-sided  derivative 
approximations  for  Px-j»  Pxg*  uxi  *  and  UX2  ™  the  directions  suggested 
by  the  corresponding  coefficients  X-j  and  X2,  the  numerical  algorithm 
will  more  accurately  model  the  physical  domain  of  dependence.  Thus, 
for  supersonic  flow,  only  upwind  information  will  enter  the  spatial 
derivative  approximations,  and  for  subsonic  flow,  a  weighted  averaging 
of  upwind  and  downwind  information  will  be  employed. 
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The  ability  of  this  scheme  to  approximate  discontinuities 
arises  from  its  accurate  modeling  of  the  physical  domain  of  dependence 
of  each  point  being  computed.  In  particular,  points  i  immediately  up¬ 
stream  of  a  shock,  where  the  flow  must  be  supersonic,  will  then  have  no 
dependence  on  the  flow  downstream  of  the  shock.  This  feature  prevents 
any  inaccuracies  that  arise  at  the  shock  from  propagating  upstream  and 
affecting  the  entire  calculation. 

The  fundamental  limitation  of  the  X-scheme  is  that  there  i- 
no  mechanism  for  an  entropy  increase  across  an  embedded  shock  since  the 
shock  is  not  treated  explicitly.  Furthermore,  the  X-scheme  is  applied 
only  to  the  continuity  and  momentum  equations:  the  energy  equation  still 
expresses  entropy  conservation  along  streamlines.  At  best,  then,  the 
X-scheme  can  be  expected  to  approximate  an  embedded  shock  as  an  isentropic 
compressive  discontinuity. 

Extending  the  X-scheme  concept  to  two-dimensional  unsteady 
flows  using  the  new  coordinate  system  requires  the  determination  oi 
characteristic  slopes  in  both  the  Z-T  and  Y-T  reference  planes.  These 


characteristic  slopes  may  be  written  as 

XZl  =  A  -  aGZ^O  +  Z^/Z^)1'2  (4.24) 

Xz2  *  A  +  aGZn(l  +  Z^/Z^2)1'2  (4.25) 

Xy-j  *  G(u -  a)  (4.26) 

Xy2  s  G(u  +  a)  (4.27) 


and  the  governing  equations  become 


PT  +  I  (xZ]pZi  +  *Z2PZ2)  +  7  (^Y-|PYi  +  ty2PY2) 

+  y(Ay2uy2  -  X y-j uy-j )/2a  +  yUz2vz2  -  Xz^zJ/Zav 
+  yGZ^uz  +  yGE  ■  o  (4.28) 

uT  +  Auz  +  £  ( X y-j UYi  +  xY2uY2)  +  GpZ^Pz/p 

+  a(Xy2PY2  -  XYiPy-j)/2y  +  GvD  =  0  (4.29) 

vj  +  \  (xZivZt  +  ^Z2VZ2)  +  Guvy 

+  a(Xz2Pz2  -  XZiPZi)/2yv  -  GuD  »  0  (4.30) 

where  v  =  (1  +  Z^2 /Z^2)1^ 2 

A  =  Z  +  GuZr  +  GvZ 
t  §  n 

D  =  v4>^  +  u<t>2 

E  =  -ud>1  +  v4>2  +  V/Gy 

The  derivative  approximations  used  in  the  X-scheme  are 


described  in  Section  5.7. 


4.3  COMPARISON  OF  AXISYMMETRIC  SHOCK-CAPTURING  PROCEDURES 

In  order  to  compare  the  relative  merits  of  the  conservation 
and  X-scheme  shock-capturing  algorithms  for  the  computation  of  flows 
over  indented  nosetips,  two  'xi symmetric  codes  have  been  developed  from 
the  analyses  presented  in  the  preceding  two  sections.  Primarily,  the 
two  codes  differ  in  their  treatment  of  the  field  points;  treatment  of 
body  and  bow  shock  points  is  essentially  the  same  for  the  two  approaches, 
and  is  equivalent  to  the  procedures  described  in  Sections  5.3  and  5.4, 

The  first  comparison  made  is  for  an  axi symmetric  indented 
shape,  the  Very  Mildly  Indented  Body  (VMIB),  on  which  wind  tunnel  tests 
have  been  conducted,  as  described  by  Reeves,  Todisco,  Lin,  and  Pall  one20. 
(Pecalibraticn  of  the  wind  tunnel  subsequent  to  the  VMIB  tests  revealed 
that  the  nominal  Mach  number  was  7.2,  and  not  8,  as  stated  in  Reference 
20.)  Figure  4.4  presents  predictions  of  bow  shock  shape  compared  to 
experimental  data  for  the  VMIB  using  the  two  schemes;  both  calculations 
used  a  16  x  28  mesh  (16  points  across  the  shock  layer  and  28  points  along 
the  body)  and  were  run  for  1000  time  steps;  convergence  criteria  (see 
Section  6.2)  were  satisfied  by  both  calculations.  This  comparison  is 
continued  in  Figure  4.5,  showing  predictions  of  the  surface  pressure 
distribution  and  the  experimental  data.  As  evidenced  by  these  figures, 
the  predictions  obtained  with  both  techniques  exhibit  good  agreement 
with  the  data.  Differences  are  apparent  in  the  two  predictions  in  the 
bow  shock  position  towards  the  downstream  boundary  of  the  noset'p  shock 
layer,  and  in  the  local  extrema  of  the  surface  pressure  distribution. 


SURFACE  PRESSURE  p/p, 


Figure  4.5.  Surface  Pressure  Predictions 
for  VMIB  at  *  7,2 
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It  should  be  noted  that,  even  though  the  VMIB  does  not  have 
a  strong  embedded  shock,  a  small  amount  of  numerical  damping  (8  *  0.008, 
with  AZ  =0.067  and  AY  =  0.085)  is  required  to  obtain  a  solution  for  this 
shape  using  the  conservation  form  of  shock-capturing.  The  differences 
noted  above  between  the  conservation  and  X-differencing  solutions  are 
entirely  consistent  with  the  effective  viscosity  implied  by  the  use  of 
numerical  damping. 

The  hazards  of  the  use  of  numerical  damping  are  illustrated 
in  Figure  4.6,  which  presents  predictions  of  the  bow  shock  shape  for  a 
sphere  in  ideal  gas  (y  =  1.4)  at  =  10  obtained  with  the  conservation 
code  using  three  values  of  the  damping  parameter  3,  with  AZ  =  0.2  and 
AY  =  0.11.  The  solution  with  6  =0  (no  damping)  agrees  well  with 
predictions  made  using  a  non-conservation  axisymmetric  code  and  the 
X-differencing  code.  The  effect  of  increased  damping  is  clearly  evident 
in  this  figure;  the  solution  with  6  *  0.02  is  fairly  close  to  that 
obtained  with  no  damping.  At  $  *  0.06,  however,  the  potential  of 
numerical  damping  to  distort  the  shock  layer  is  graphically  illustrated. 

For  most  calculations,  of  course,  a  comparison  such  as  that 
in  Figure  4.6  is  not  possible,  and  the  risk  of  obtaining  a  "reasonable" 
conservation  solution  with  numerical  damping  that  is  a  poor  approximation 
to  the  inviscid  flow  being  computed  is  great.  In  this  regard,  the 
X-differencing  approach  is  judged  to  be  superior  to  the  conservation 
approach  in  that  the  X-differencing  scheme  eliminates  the  uncertainties 
inherent  in  the  numerical  damping  required  to  generate  solutions  with  the 
conservation  formulation. 
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Figure  4. 


^  A 


.  Effect  of  Damping  on 
Shock  Shape  Predictions 
with  Conservation  Scheme 
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Another  comparison  of  the  two  shock-capturing  approaches 
under  consideration  has  been  made  for  the  Mildly  Indented  Body  (MIB), 
another  of  the  axisynmetric  shapes  considered  by  Reeves,  Todisco,  Lin, 
and  Pal lone20.  Figure  4.7  shows  the  bow  shock  shape  for  this  indented 
body  (on  which  an  embedded  shock  does  form)  as  predicted  by  the  two 
shock-capturing  schemes  compared  to  the  data  obtained  at  ^  *  9.  (In 
Reference  20,  the  Mach  number  for  the  MIB  tests  was  given  as  12; 
recalibration  of  AEDC  Tunnel  F  has  indicated  the  true  Mach  number  was 
approximately  9.)  The  X-scheme  has  produced  a  converged  solution  for 
this  case  (1000  steps),  while  the  conservation  solution  shown  (400  steps) 
has  not  converged;  in  fact  the  oscillations  arising  in  the  conservation 
solution  from  the  smeared  embedded  shock  are  growing  and  the  solution 
is  diverging.  Control  of  these  oscillations  requires  application  of  a 
level  of  numerical  damping  that  grossly  distorts  the  predicted  bow 
shock  shape. 

As  is  evident  from  Figure  4.7,  the  X-differencing  scheme 
produces  fair  agreement  with  the  data,  but  locates  the  triple  point 
(where  the  embedded  shock  intersects  the  bow  shock)  too  far  downstream 
and  the  downstream  bow  shock  too  far  away  from  the  body.  These 
discrepancies  arise  because  of  the  lack  of  accurate  modeling  of  the 
shock  Intersection  point  and  the  assumption  of  isentropic  flow  across 
the  embedded  shock. 

In  an  attempt  to  combine  the  beneficial  features  of  both  the 
X-differencing  and  conservation  approaches  to  embedded  shock  predictions, 
conservation  calculations  were  made  for  the  MIB  using  second-order 
accurate  upwind  differences  in  regions  of  supersonic  flow.  Although 
slightly  reducing  the  amount  of  numerical  damping  required  to  obtain  a 
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Figure  4,7,  Shock  Shape  Prediction*  for  fllB  at  K,  = 


solution,  this  modified  conservation  approach  did  not  produce  any 
significant  improvement,  in  the  conservation  results  depicted  in  Figure 
4.7. 

This  case  serves  as  an  example  of  the  ability  of  the 
X-differencing  scheme  to  compute  inviscid  time-dependent  flows  with 
embedded  shocks  in  a  stra»ght-forward  manner  without  resorting  to  the 
artifice  of  numerical  damping,  thus  avoiding  the  problem  of  inconsistency 
of  the  damped  finite  diff  ance  equations  with  the  inviscid  flow 
equations. 

It  should  be  noted  that  frequently  many  attempts  are  required 
to  produce  a  “gr  i"  conservation  solution,  the  problem  being  to  find  the 
smallest  possible  amount  of  mping  that  will  sufficiently  control  the 
oscillations  to  allow  the  solution  to  proceed  without  failure,  without 
overly  distorting  the  flow  being  computed.  In  some  cases,  no  such  value 
of  the  damping  parameter  has  been  found  for  the  conservation  procedure 
developed  in  this  research. 

Based  on  a  number  of  comparisons  of  the  X-differencing  and 
conservation  solutions,  of  which  the  comparisons  presented  in  this 
section  are  cniy  a  sample,  it  is  concluded  chat  the  X-differencing  scheme 
offers  a  significant  adW.tage  over  the  conservation  solution  of  the 
ambeuded  shock  problem.  Since  the  X-scheme  not  produce  oscillations 
in  the  solution  as  it  "captures"  on  embedded  sho^;-,  as  does  the  con¬ 
servation  formulation,  numerical  damping  is  not  required,  and  the  question 
of  consistency  of  the  prob.em  being  solved  numerically  with  the  inviscid 
problem  does  not  arise.  Furthermore,  the  X -scheme  has  been  found  to  oe 
far  more  efficient  to  apply  to  a  g  ■  ven  problem,  general  ly  requiring  only 
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one  calculation,  where  many  conservation  solution  attempts  are  typically 
required  to  determine  the  optimum  value  of  the  damping  parameter.  Be¬ 
cause  of  these  comparisons,  the  X-differencing  scheme  is  the  approach 
to  the  embedded  shock  problem  that  has  been  selected  for  extension  to 
the  three-dimensional  problem. 

It  must  be  remembered,  however,  that  the  X-differencing 
scheme  can,  at  best,  only  approximate  the  structure  of  flows  with  em¬ 
bedded  shocks  because  of  its  assumption  of  isentropic  flow.  As  the 
strength  of  the  embedded  shock  grows,  this  approximate  embedded  shock 
solution  can  be  expected  to  become  less  accurate.  Improving  the  aco'racy 
of  embedded  shock  calculations  will  require  strict  enforcement  of  the 
Rankine-Hugoniot  conditions,  necessitating  a  shock-fitting,  instead  of 
a  shock-capturing,  approach. 

4.4  X-DIFFERENCING  SCHEME  IN  THREE  DIMENSIONS 

The  extension  of  the  axi symmetric  X-differencing  scheme 
described  in  Section  4.2  to  three  dimensions  requires  the  determination 
of  characteristic  slopes  in  the  Z-T,  Y-T,  and  X-T  reference  planes. 

These  six  characteristic  slopes  may  be  written  as: 


Xz-j  =  A  -  aGZ^v^ 

(4.31) 

Xz«  =  A  +  aGZ^v^ 

(4.32) 

Xy-j  *  B  -  aGY^v^ 

(4.33) 

XY2  *  B  +  aGY^Vy 

(4.34) 

XXl  =  XQ(w-a)/y 

(4.35) 

XX2  s  X0(w+a)/y 

(4.36) 

88. 


where 


-I’"*  V'V  +  <Vc/Zn  +  "♦  +  VV^V)''1 

vy  *  0  +  (5^  ♦  Y0/Vc)!/6Vl,,J 

In  the  X-differencing  scheme  the  continuity  and  momentum 
equations  described  in  Section  3.4  (Equations  ( 3.69 )- ( 3 . 72 ) )  are  re 
written  in  a  form  that  allows  accurate  modeling  of  the  domain  of 
dependence  of  each  point.  The  resulting  three-dimensional  time- 
dependent  equations  become 

PT  +  1  <xZlPZl  +  XZ2PZ2>  +  J  (xYipYi  +  XY2PY2) 

+  \  (xXiPXi  +  XX2PX2)  +  7Y(Xy2>JY2  "  xYiuYi)/avY 
+  ^y(Xz2vZ2  ‘  xZi  VZ]  )/avz  +  7  y(xX2wX2  ‘  xX]wXi  )/a 
+  yG(Z^uz  +  E  +  V/Gy)  +  y!(^Y?  +  YQ)w  y 

+  (Ve  +  Vn  +  Ze)wZ]/y  s  0 

u-j-  +  Auz  +  -j-  (?. y-j  Uyi  +  Xy^Jyg)  +  ^ux  + 

+  vwF/y  -  $w2/y  +  GpZ^/p 
+  \  a  (Xy2pY2  "  xYipYi)/YVy  *  0 

VT  +  1  (xZ]vZl  +  XZ2VZ2)  +  Bvy  +  Cvx  -  uGD 

-uwF/y  -  Cw2/y  +  a  (Xz2Pz2  “  xZipZi)/wz  a  0 
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(4.37) 


(4.38) 


(4.39) 


(4.40) 


wT  +  Awz  +  Bwy  +  j  (XxiwXi  *'  XX2WX2^  +  w(u$  +  v^/y 

+  P  +  V*  +  Z0)PZ  +  (Vc  +  VPYJ/py 
+  \  a  (XX2PX2  *  xX1pXt)/y  *  0 

The  coefficients  appearing  in  the  above  equations  are  those  defined 
in  Section  3.4. 

As  with  the  axisynmetric  version  of  the  X-scheme,  the  energy 
equation  (conservation  of  entropy  along  streamlines)  is  unaltered  in 
the  X-differencing  formulation.  Centerline  points  are  treated  using 
the  special  forms  of  the  governing  equations  derived  in  Section  3.6. 
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SECTION  5 


NUMERICAL  PROCEDURES 

5.1  TIME-ASYMPTOTIC  SOLUTION  &  INITIALIZATION 

The  approach  selected  for  the  solution  of  this  steady  inviscid 
flow  problem  is  the  time-asymptotic  relaxation  method.  In  this  approach, 
the  steady  flow  solution  is  sought  as  the  limit  of  a  time-dependent  flow 
with  time  invariant  boundary  conditions  as  t  -*■  «.  Because  the  time- 
dependent  equations  are  hyperbolic,  this  method  allows  solution  of  the 
steady  flow  problem,  which  is  a  boundary  value  problem,  as  a  mixed 
initial -boundary  value  problem  with  a  forward-marching  (in  time) 
numerical  technique. 

Implicit  in  this  approach  is  the  assumption  that  a  steady 
flow  limit  exists  and  is  unique.  Purely  unsteady  flows,  with  an 
oscillating  bow  shock,  can  occur  on  severely  indented  nosetips,  as 
analyzed  by  Reeves28.  In  principle,  a  time-dependent  numerical  technique 
can  compute  such  oscillatory  flows.  However,  such  flows  are  dominated 
by  viscous  separation  effects,  and  thus  are  not  considered  in  the  in¬ 
viscid  problem  being  addressed. 

To  properly  pose  the  initial  value  problem,  it  is  necessary 
to  completely  define  the  flow  field  at  some  arbitrary  instant  of  time, 
t  ■  0.  Theoretically,  any  estimate  of  the  flow  field  at  this  instant 
will  suffice,  with  the  solution  eventually  converging  to  the  (assumed) 
unique  steady  state.  In  practice,  however,  it  is  best  to  use  as 
reasonable  an  initial  flow  field  estimate  as  possible,  for  grossly 
inaccurate  initializations  produce  large  gradients  in  the  flow  that  the 
numerical  scheme  is  incapable  of  handling.  Additionally,  estimates  of 
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the  Initial  state  that  are  reasonable  estimates  of  the  final  steady 
flow  field  tend  to  accelerate  convergence  of  the  time-dependent  technique. 

Initial  data  is  generated  for  this  nosetip  flow  field  technique 
by  neglecting  the  effects  of  circumferential  flow  variation.  First,  an 
estimate  of  the  shock  stand-off  distance  at  the  centerline  Is  obtained 
from 


Aq  «  Rn  [0.61 37/ (Mj-1)  +  0.131  (5.1) 

which  is  a  curve  fit  of  numerical  results  for  spheres  at  zero  angle  of 
attack.  From  this  shock  point  at  the  centerline,  estimates  of  the  shock 
shape  in  each  meridiona.l  plane  can  be  obtained  using  a  correlation  for 
shock  slope  on  a  sphere  in  terms  of  the  equivalent  body  angle  9b  from 
Abbett  and  Davis29 as: 

a  *  0.5236  +  0.33339b  +  0.2122eb2  (5.2) 

where  the  equivalent  body  angle  is  defined  as 

3yfa 

eb  ■  tan-1  [  -  a,  cos<}>  -  B  sin<J>  ]  (5.3) 

thus  including  some  effects  of  angle  of  attack  and  sideslip.  Strictly 
speaking,  the  correlation  given  by  Equation  (5.2)  applies  along  a  body- 
normal;  for  this  application  it  is  assumed  valid  along  K  =  constant 
curves. 

Given  all  of  the  shock  points,  tiie  shock  slopes  c^  and  cg 
may  then  be  evaluated  using  finite-difference  expressions.  From  the 
shock  slopes  and  f^eestream  conditions  the  downstream  shock  properties 
are  then  evaluated  from  the  Rankine-Hugoniot  conditions. 


Next,  the  flow  conditions  at  the  body  are  approximated.  From 
the  modified  Newtonian  impact  theory,  the  surface  pressure  may  be  ex¬ 
pressed  as 


p  =  pQ  sin2efa  (5.4) 

where  the  equivalent  body  angle  is  given  by  Equation  (5.3)  and  pQ  is  the 
normal  shock  stagnation  pressure.  Assuming  that  the  normal  shock  value 
of  entropy  applies  at  the  body  surface  (see  Section  3.1),  the  static 
enthalpy  h(p,s)  can  then  be  evaluated  from  the  thermodynamic  state 
relations.  Since  the  total  enthalpy  H  is  constant  for  a  steady  inviscid 
flow  and 

H  =  h  +  ^  q2  -  h  +  £  (u2  +  v2  +  w2)  (5.5) 

the  total  velocity  V  at  the  body  surface  can  be  found.  Neglecting 
crossflow  (w  =  0)  and  imposing  the  kinematic  boundary  condition  results 
in 

u  =  V/(l  +  b?2)l/z  (5.6) 

v  =  ub^  .  (5.7) 

To  complete  the  specification  of  the  initial  flow  field, 
linear  distributions  of  P,s,v,  and  w  in  n  are  assumed  between  the  body 
and  shock  along  £  =  constant  lines.  The  final  velocity  component  u  is 
assigned  from  the  total  enthalpy  relation.  Equation  (5.5). 
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From  this  assumed  initial  state  the  time-dependent  equations 
can  be  integrated  forward  in  time  to  the  steady  state  solution.  The 
problem  still  remains,  however,  of  determining  when  a  solution  is 
sufficiently  converged  to  be  considered  an  adequate  representation  of 
the  steady  solution. 

Convergence  of  an  iterative  solution  can  be  defined  in  terms 
of  variations  of  the  dependent  variables;  i.e.,  when  the  percentage 
change  from  step  to  step  of  all  dependent  variables  at  all  mesh  points 
is  less  than  some  e,  the  solution  can  be  assumed  converged.  The  nature 
of  the  inviscid  nosetip  time-dependent  problem  does,  however,  con¬ 
veniently  provide  other  parameters  that  can  be  used  as  a  measure  of 
convergence. 

During  the  course  of  a  time-dependent  calculation  the  bow 
shock  position  continuously  adjusts  from  its  assumed  initial  position 
to  its  final  (converged)  steady  state  position.  Thus,  the  root-mean- 
square  of  the  velocities  of  the  shock  points  used  in  the  calculation 
serves  as  a  convenient  measure  of  the  convergence  of  the  solution. 
Additionally,  the  motion  of  the  bow  shock  during  the  calculation  provides 
another  criterion  for  convergence:  the  conservation  of  total  enthalpy. 
Unlike  the  steady  case,  the  total  enthalpy  is  not  conserved  across  a 
moving  shock;  thus  the  total  enthalpy  will  vary  throughout  the  flow 
field  during  a  time-dependent  calculation.  The  variation  of  the  total 
enthalpy  from  its  known  freestream  value  then  serves  as  a  direct 
measure  of  convergence  of  the  calculation,  since  this  difference  in  total 
enthalpy  will  diminish  only  if  all  shock  velocities  are  diminishing. 
Details  on  the  actual  convergence  criteria  used  to  judge  the  merits  of 
a  calculation  are  provided  in  Section  6.2. 
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5.2  NON-CONSERVATION  FIELD  POINT  TREATMENT 

At  field  points  (i.e.,  all  computational  grid  points  not 
located  either  on  the  body  or  the  shock),  the  governing  time-dependent 
equations  in  non-conservation  form  are  solved  using  the  MacCormack  7 
finite-difference  scheme,  with  some  variations.  This  widely  used 
numerical  technique  is  a  two-step  predictor-corrector  scheme  and  main¬ 
tains  second  order  accuracy  in  both  mesh  spacing  and  time. 

To  illustrate  the  MacCormack  technique,  consider  the  simple 
hyperbolic  equation 

ft  +  Afx  =  0  (5.8) 

k 

where  f  may  represent  either  a  scalar  or  a  vector  quantity,  and  fn 
represents  f(tQ  +  kAt,  xQ  +  nAx).  In  the  predictor  stage,  the  spatial 
derivative  is  approximated  as  a  forward  difference 

fx  ■  <f£+i  -  (5-9> 

~k 

and  the  predicted  value  of  fn  at  t  =  t  +  (k+l)At,  denoted  by  fn,  is 
obtained  from 

■  fk  -  A  4t(fk+,  -  f*)Mx  .  (5.10) 

In  the  corrector  stage,  f  is  approximated  using  backward 
differences  of  the  predicted  data: 

fx  •  -  *£,>/*«  •  (5-n) 
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Overall  second  order  accuracy  is  achieved  by  obtaining  the 
final  corrected  value  from 

fnk+1  -  \  [fk  +  fk  -  A  At(fk  -  f^J/Ax]  .  (5.12) 

Assuming  the  coefficient  (or  coefficient  matrix)  to  be  constant,  a 
truncation  error  analysis  shows  the  leading  order  error  of  this  scheme 
to  be  -1/C  [fttt(At)2  +  Afxxx(Ax)2]  for  the  simple  equation  considered 
here. 

The  standard  MacCormack  scheme  is  not  used  consistently 
throughout  non-conservation  calculations,  however.  As  noted  by  Moretti17, 
the  standard  scheme  is  inappropriate  for  the  approximation  of  convective 
derivatives;  i.e.,  Lagrangian  derivatives  should  be  approximated  using 
upwind  differences  only.  To  maintain  second  order  accuracy  for  these 
cases  in  the  context  of  a  predictor-corrector  scheme,  the  standard 
MacCormack  scheme  can  be  modified  by  replacing  Equation  (5.9)  with 

fx  “  -  3fn-1  +  fn-2>/4*  (5J3) 

when  backward  differences  are  required,  or  by  replacing  Equation  (5.11) 
with 

fx  ■  <-t2  +  3t,  -  **>«*  (5-14> 

when  forward  differences  are  required.  Use  of  either  of  the  above 
modifications  retains  the  overall  order  of  accuracy  of  the  calculation 
with  a  truncation  error  of  1/6  fttt(At)2  +  V3  ^xxx^x^2, 
convective  differencing  version  of  the  MacCormack  predictor-corrector 
scheme  is  used  in  the  non-conservation  form  of  the  governing  equations 
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for  the  entropy  derivatives  In  Equation  (3.73)  and  for  the  derivatives 
in  the  convective  terms  in  the  momentum  Equations  (3.7Q)-(3.72). 

Because  the  MacCormack  scheme  is  explicit,  limitations  on 
the  time  step  to  be  taken  must  be  Imposed  to  prevent  unstable  calcula¬ 
tions.  According  to  the  Courant-Friedrichs-Lewy  (CFL)30criterion  for 
stability,  the  time  step  must  be  selected  such  that  the  numerical  domain 
of  dependence  does  not  exceed  the  physical  domain  of  dependence  of  the 
point  being  computed.  The  allowable  time  step  may  be  evaluated  as  the 
minimum  value  at  all  points  of  At  computed  as 

At  *  min(Ax,Ay)/ [/2{(u2  +  v2  +  w2)1/,2+  a)]  .  (5.15) 

This  form  assumes  that  the  circumferential  spacing  (yA<f> )  will  not  be  the 
controlling  length  scale,  a  reasonable  assumption  for  the  ablated  nosetip 
geometries  usually  encountered. 

5.3  BODY  POINT  TREATMENT 

Special  computational  procedures  are  required  at  the  boundary 
points  of  the  computational  region  (i.e.,  body  points  and  bow  shock 
points),  where  the  governing  partial  differential  equations  must  be 
solved  in  conjunction  with  the  boundary  conditions  that  are  to  he  imposed. 

At  the  body  points,  the  Kentzer-Moretti  predictor-corrector 
scheme  is  used  to  accurately  model  the  physics  of  the  flow  at  the 
impermeable  boundary.  In  this  scheme,  the  discretization  of  the  boundary 
conditions  suggested  by  Kentzer”  was  extended  to  the  predictor-corrector 
format  by  Moretti  and  Pandolfi9  .  The  parallels  between  this  scheme  and 
the  numerical  method  of  characteristics  solution  at  the  body  point  have 
been  examined  by  Hall14. 
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In  the  application  of  the  Kentzer-Moretti  scheme  to  the 
three-dimensional  time-dependent  problem,  the  continuity  relation  given 
by  Equation  (3.69)  is  replaced  by  the  characteristic  compatibility 
condition  expressed  by  Equation  (3.99).  The  circumferential  momentum 
relation.  Equation  (3.72),  is  used  directly,  and  the  surface  entropy  is 
assigned  as  the  normal  shock  value  of  the  entropy,  as  discussed  1r. 

Section  3.1.  (DeNeef31  has  recently  published  an  efficient  method  of 
coding  the  compatibility  condition.) 

To  complete  the  specification  of  the  body  point  procedure, 
two  more  relations  are  required. 

The  kinematic  boundary  condition.  Equation  (3.97),  provides 
one  of  these  relations.  The  other  may  be  obtained  from  a  linear  combina¬ 
tion  of  the  E-  and  n-  momentum  equations,  given  in  computational  coordi¬ 
nates  by  Equations  (3.70)  and  (3.71).  Defining  a  "tangential"  velocity 
component  v  as 

v  =  u  +  b^  v  (5.16) 

a  time-dependent  "tangential"  momentum  equation  may  be  written  as 

vT  -  uT  +  b^  vT  (5.17) 

with  Uj  and  Vy  being  evaluated  from  Equations  (3.70)  and  (3.71). 

The  solution  at  a  body  point  is  obtained  by  solving  Equations 
(3.99),  (3.72),  and  (5.17)  with  the  standard  MacCormack  scheme,  except 
that  Z-derivatives  must  always  be  approximated  as  forward  differences, 
requiring  the  use  of  Equation  (5.14)  in  the  corrector  stage.  Once  new 
values  of  v  and  w  have  been  computed,  the  other  velocity  components  can 


98. 


be  determined  from  Equation  (5.16)  and  the  kinematic  boundary  condition 
to  be 


v  n 


(5.18) 


u  «  v  -  vb 


5 


(5.19) 


It  was  discovered  during  the  formulation  of  this  body  point 
procedure  that  the  form  of  the  "tangential"  momentum  equation  described 
above,  where  vT  is  a  linear  combination  of  Uj  and  vT  determined  from 
Equations  (3.70)  and  (3.71),  respectively,  must  be  used  when  convective 
Y-differences  are  taken.  In  Moretti's17  original  axisymnetric  formula¬ 
tion  of  this  scheme.  Equation  (5.17)  was  written  as 

Vj  *  B(vy  -  vb^)  +  .  .  .  (5.20) 


which  follows  analytically  from 

vvv  vb« 


(5-21) 


In  theory,  this  approach  is  entirely  equivalent  to  tha  result  obtained 
from  the  linear  combination  of  Equations  (3.70)  and  (3.71),  namely 

vT  *  Buy  +  b^BVy  +  .  .  .  .  (5.22) 

However,  when  convective  Y-differences  are  used  for  Vy  in  Equation  (5.20) 
and  for  uY  and  Vy  in  Equation  (5.22),  the  resulting  finite  difference 
expressions  are  not  equivalent.  A  truncation  error  analysis  shows  that 
second  order  accurate  convective  differencing  for  Vy  at  Y  *  mAY  produces 


Vy  ■  Uy  + 


(5.23) 


vv  +  vbrr  +  O(AY) 
m-1  Y 

which  does  not  produce  an  accurate  representation  of  Uy  +  b^Vy.  Thus, 
the  form  for  Vy  given  by  Equation  (5.17)  is  the  preferred  form  for  this 
approach,  and  has  been  shown  to  produce  more  accurate  solutions,  in  the 
sense  of  producing  smaller  errors  in  total  enthalpy  at  the  body  when 
the  solution  has  converged  to  the  steady  state. 

5.4  BOW  SHOCK  POINT  TREATMENT 

The  Kentzer-Moretti  predictor-corrector  scheme,  described  in 
the  previous  section,  Is  also  used  for  the  computation  of  bow  shock 
points.  Unlike  the  body  point  procedure,  however,  the  use  of  this  scheme 
at  bow  shock  points  produces  an  equation  for  the  shock  acceleration 
derived  from  a  characteristic  compatibility  condition,  which  can  be 
integrated  twice  in  time  to  compute  the  shock  velocity  and  position. 

The  procedure  outlined  below  is  essentially  that  developed  and  used 
extensively  by  Moretti ;  e.g.,  as  in  Reference  17;  application  of  this 
method  has  recently  been  simplified  by  deNeef31. 

Derivation  of  the  equation  for  the  shock  acceleration  starts 
with  the  characteristic  compatibility  condition  given  by  Equation  (3.95). 
The  normal  frees tream  velocity  component  relative  to  the  moving  shock 
may  be  written  as 

5.  -  u00N1  +  (v.  -  cT/G)N?  +  w08N3  (5.24) 
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where 


u  -  V  [C  cosBcosa  +  S(s1nBs1n6  +  cosBsInocosB)) 

SB  • 


(5.25) 


v  •  V  I-S  cosBcosa  +  C(s1nBs1n6  +  cosB$1nacose)l 


(5.26) 


w„  ■  VfsinBcose  -  cosBsinasIne) 


(5.27) 


and  a  and  B  are  the  angles  of  attack  a^  sideslip,  respectively.  The 
unit  normal  to  the  shock  Is 


ns  «  N1  i  +  N2  j  +  N3  k 


(5.28) 


where 


N.j  *-6c^/v 


N2  *  G/v 


(5.29) 


(5.30) 


N3  “  (^  -  V5'C9)/yV 


(5.31) 


v  «  (G2(l  +  cc2)  *  -  cQ)2/y2] 1/2 


(5.32) 


The  normal  downstream  velocity  component  relative  to  the  moving  shock 


u  *  uN-j  +  (v  -  Cy/GjNg  +  wNj 


and  the  downstream  velocity  components  may  be  written  as 


u  •  u„  +  (u  -  u#0)N1 

(5.34) 

v  ■  v.  ♦  (u  -  uJN2 

(5.35) 

w  -  ww  +  (u  -  uJN3 

(5.36) 

The  time  derivatives  that  appear  In  Equation  (3.95)  can  now  be  replaced 
with  expressions  Involving  the  shock  acceleration  by  differentiating  the 
Rankine-Hugoniot  conditions  to  obtain  PT,  iLp  and  u-p  and  thus  also 
u-p,  vT,  and  wT  from  Equations  (5.34)-(5.36) . 

The  final  equation  for  the  shock  acceleration  may  then  be 

written  as 

4 

cn  3  CD  tiRi  -  T2  -  T3]/T1  (5.37) 

t*l 

where 

T^  *  Cj  +  i^Cg  ^  ^3^11  ^4^1 3  (5.38) 

T2  3  X[PZ  +  i-2uz  +  i-3vz  +  t4wz)  (5.39) 

T3  =  C4  +  S'2C10  +  *3C12  +  £4C14  *  ^5’40^ 

The  Z-derivatives  in  Equation  (5.39)  must  be  evaluated  with  backward 
differences,  using  Equation  (5.13)  in  the  predictor  calculation.  The 
coefficients  in  the  above  expressions  are  derived  in  the  Appendix. 


•  / 
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In  the  predictor  stage ,  Equation  (5.37)  can  be  Integrated 

to  yield 

cT  *  Cj  ♦  Cjj  AT  (5.41) 

c  *  ck  +  c*  AT  .  (5.42) 

The  final  values  of  the  shock  speed  and  position  are  obtained  in  the 
corrector  stage  of  the  calculation  from 


:k+^  “  7  (ck  +  Cj  +  C-q-  &T) 

(5.43) 

:k+1  -  \  (Ck  ♦  C  ♦  Ck+1  AT) 

(5.44) 

Once  the  new  •‘hock  positions  have  been  determined  (predictor 
or  corrector),  the  shock  slopes  c^  and  cQ  may  be  evaluated  from 

c^  *  Cy  (5.45) 

and 

c0  =  CX  X0  +  CY  Y9  (5.46) 

where  cY  and  cx  are  evaluated  from  difference  formulas.  At  most  shock 
points,  Cy  is  evaluated  at  Y  =  mAY  from  a  standard  four  point  formula  as 

CY  *  (cm-2  '  6cm-l  +  3cm  +  2cm+l)/6AY 


(5.47) 


At  the  last  two  points  In  each  6-plane»  the  expression 

CY  *  <-2cm-3  +  9S-2  ’  18cm-l  +  llcm>/siY  <5‘48) 

is  used,  while  at  the  point  adjacent  to  the  centerline  the  form  used  is 

CY  ■  <-Vl  -3cm  +  6Vl  -  W'6AY  <5'49> 

c^  is  evaluated  at  X  ■  AAX  using  the  centered  two-poir.t  difference 
formula: 


Cx  *  ^j,+l  ”  c^_i)/2AX  .  (5.50) 

The  other  shock  derivatives  required,  c^T  and  cQT,  are 
determined  using  the  difference  formulas  given  above,  but  using  cT,  the 
shock  velocity,  in  place  of  c. 

Knowing  the  shock  slopes  and  velocities,  the  properties  down¬ 
stream  of  the  shock  are  determined  directly  from  the  Rankine-Hugoniot 
conditions: 

pu  *  Pj5»  (5.51) 


pu2  =  p.  +  pjj.2 

(5.52) 

u2  =  1  V 

(5.53) 

For  an  ideal  gas,  these  equations  can  be  solved  analytically  for  p,  p, 
and  u;  for  an  equilibrium  real  gas  an  iterative  solution  is  required. 
To  avoid  this  time-consuming  Iteration  for  real  gas  calculations,  a 
table  of  shock  properties  as  a  function  of  u,,,  is  created  at  the  start 
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of  the  calculation,  thus  requiring  only  interpolation  to  determine  the 
downstream  properties  at  any  shock  point  during  the  computation.  Once 
u  has  been  determined,  the  downstream  velocity  components  can  be  computed 
from  Equations  (5.34)-(5.36) . 

5.5  CENTERLINE  POINT  TREATMENT 

Modifications  of  the  numerical  procedures  described  in  the 
preceding  sections  are  required  for  computational  grid  points  located 
on  the  singular  centerline.  These  special  procedures  required  at  the 
centerline  are  made  necessary  by  the  special  forms  of  the  governing 
equations  derived  in  Section  3.6,  in  which  the  equations  at  the  center- 
line  are  written  with  terms  evaluated  only  in  the  planes  <J>  =  0 ,  tt/2,  it, 
and  3tt/2. 

Within  any  4>  plane,  only  forward  differences  in  Y  are  possible 
at  the  centerline.  In  order  to  maintain  a  two-sided  predictor-corrector 
sequence  for  the  Y-differences  at  the  centerline,  terms  in  the  governing 
equations  are  evaluated  in  the  <t>  =  0  and  $  =  tt/2  planes  in  the  predictor 
stage,  and  in  the  4>  *  tt  and  <J>  =  3tt/2  planes  in  the  corrector  stage.  The 
resulting  numerical  scheme  thus  utilizes  information  from  all  directions, 
in  the  spirit  of  the  finite  difference  scheme  used  at  points  not  located 
on  the  centerline. 

At  field  points  along  the  centerline,  the  equations  to  be 
solved  with  this  procedure  are  given  by  Equations  (3.1 08 )-(3.112) .  A 
difficulty  arises,  however,  in  the  application  of  this  technique  when 
the  coordinate  transformation  to  (£,n,e)  space  is  not  axisynmetric.  In 
this  case,  variations  in  the  complex  scaling  factor  g  along  the  center- 
line  among  the  <j>  planes  result  in  the  non-correspondence  of  the 


105. 


computational  grid  points  in  physical  space.  In  other  words,  even 
though  the  computational  grid  points  at  the  centerline  are  equally 
spaced  in  Z  for  all  4>  planes,  grid  points  with  the  same  Z  value  in 
different  <J>  planes  may  represent  different  points  in  physical  space. 

These  differences  are  minimized  by  the  use  of  the  stretching  parameter 
a^  in  the  coordinate  transformation  defined  in  Section  3.2,  and  thus 
linear  interpolation  procedures  can  be  used  to  resolve  the  differences 
without  any  noticeable  effect  on  accuracy. 

Calculations  of  the  body  and  shock  points  along  the  center- 
line  do  not,  of  course,  experience  this  difficulty,  thus  simplifying 
the  numerical  approach  at  these  points.  Details  on  the  modifications 
required  for  the  characteristic  slopes  and  compatibility  conditions  at 
these  boundary  points  are  provided  in  Section  3.6. 

A  special  procedure  is  required  for  the  calculation  of 
shock  slopes  at  the  centerline,  where  accuracy  requires  a  formulation 
that  utilizes  information  from  more  than  one  <j>  plane.  In  particular, 
the  calculation  of  c^  in  some  plane  also  requires  information  about 
the  shock  shape  from  the  plane  +  tt.  A  standard  difference  procedure 
cannot  be  used,  however,  when  an  asymmetric  mapping  is  invoked;  i.e., 
when  the  mesh  spacings  AC  and  the  transformed  shock  locations  c  are 
not  necessarily  equal  in  the  4>1  and  4>1  +  tt  planes. 

At  the  centerline,  S  =  1  and  C  =  0  for  all  <J>  planes.  It 
then  follows  that,  with  (  )-|  denoting  the  <}>.|  plane  and  (  )2  denoting 

the  <|>i  +  ir  plane, 

(dn/G)1  =  (dn/G)2  (5.54) 
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and,  in  the  limit  as  y+0, 


(dC/G)]  *  (d5/S)2  .  (5.55) 

Using  these  relations,  a  four- point  difference  formula  for  can  be 
derived  for  unequally  spaced  base  points; 

c^  =  Acq  +  Bc^  +  Cc2  +  Dc_i  (5.56) 

where 

Cg  =  C  (  0,(j>-|  ) 

C-|  =  c(A£^  » d> i  ) 

c2  =  c(2A£-|  j^i  ) 

c_i -  Cg  +  [c(A^2 *4>-j  )  ~  c(0,<j>.j  +  it)]  G-j /G2 

Defining 

a  =  G1A?2/G2A51  (5.57) 

where  Af^  and  A^2  are  the  mesh  spacings  in  the  respective  <p  planes,  the 
coefficients  in  Equation  (5.56)  are  given  by 

A  =  (a2-5a  +  2)/2A^a(a+l )  (5.58) 

B  =  2/AC-|  (a+1 )  (5.59) 

C  =  -a/2A?1(a+2)  (5.60) 

D  =  2/A?1a(a2  +  3a  +  2) 


(5.61) 


The  derivative  c^j  is  evaluated  at  the  centerline  using  the  same  approach. 
For  the  case  of  an  axisymmetric  mapping  with  ct  =  1,  the  difference  formula 
given  by  Equation  (5.56)  reduces  to  the  standard  four-point  formula  given 
by  Equation  (5.47). 

Once  the  solution  at  a  centerline  point  has  been  obtained 
within  any  one  $  plane,  the  solution  within  other  <f>  planes  will  also 
have  been  determined.  Thermodynamic  properties,  such  as  P  and  s,  corre¬ 
spond  Erectly,  but  a  change  in  the  value  of  <J>  requires  a  rotation  of 
the  velocity  components.  The  transformed  velocity  components  within 
any  «+>  plane  can  be  written  in  terms  of  the  known  components  in  some 
plane  4>  =  from 

u(<j>)  =  u(4>i )  [cos<)>cos<j>i  +  sin<j>sin<j>.j] 

+  w(<J>1 )  [sin4>cos4>-j  -  cos<j>sin<j>.j] 

v(4>)  =  v(<j>1) 

w(4>)  =  u(d>-j )  [cosefrsind)^  -  si n<j>cos<t>1  ] 

+  w(4>1 )  [cos^cosd)^  +  sin^sin^]  (5.64) 

Similarly,  the  shock  slope  c^  at  the  centerline  within  any  <j>  plane  may 
be  expressed  in  terms  of  the  slopes  in  the$  =  0  and<J>  -  tt/2  planes  as 

c^(<J>)  =  c^(<j>=0)cos(j>  +  cc(d>=f)sin<|>  .  (5.65) 


(5.62) 

(5.63) 
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The  transformed  centerline  shock  velocities  may  also  be  related  from 
the  condition  that  the  quantity  Cy/G,  which  corresponds  to  the  physical 
shock  velocity,  is  the  same  for  all  <j>  planes;  thus 

Cy(<M  a  G(d> )cy(<t>y  )/G(d>^ )  .  (5.66) 

5.6  CONSERVATION  FIELD  POINT  TREATMENT 

The  conservation  law  approach  to  shock-capturing,  described 
in  Section  4.1,  is  one  of  two  methods  investigated  for  the  calculation 
of  embedded  shocks  within  inviscid  shock  layers  on  reentry  vehicle 
nosetips.  This  approach  requires  the  solution  of  the  governing  equations 
written  in  conservation  form  at  all  field  points  using  a  version  of  the 
MacCormack  finite  difference  scheme  described  in  Section  5.2.  After  each 
integration  step,  the  conservation  solution  may  be  damped  and  smoothed, 
as  required. 

One  change  has  been  made  in  the  MacCoimack  scheme  for  appli¬ 
cation  to  the  conservation  calculations.  Because  the  oscillations  inherent 
in  conservation  calculations  near  a  discontinuity  arise  from  approximating 
derivatives  across  the  discontinuity,  the  modification  made  to  the 
MacCormack  scheme  is  designed  to  minimize  the  number  of  differences  taken 
across  the  discontinuity.  This  is  accomplished  by  more  closely  aligning 
the  numerical  domain  of  dependence  of  any  given  field  point  with  the 
probable  orientation  of  an  embedded  shock  within  the  computational  mesh. 

To  illustrate  this  point,  consider  the  numerical  domain  of 
dependence  of  a  point  Y  *  mAY  and  Z  =  nAZ.  Defining  the  numerical  domain 
of  dependence  of  a  point  as  all  grid  points  that  can  effect  the  solution 
at  that  point  in  one  time  step,  the  standard  MacCormack  scheme  produces 
the  domain  of  dependence  (in  the  Y-Z  plane)  shown  in  Figure  5.1. 
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Figure  5.1.  Numerical  Domain  of 

Dependence  for  Regular 
MacCormack  Scheme 

Since  the  Y-coordinate  direction  is  nearly  normal  to  the  body 
an  embedded  shock  can  be  expected  eo  be  aligned  within  the  Y-Z  mesh  as 
shown  in  Figure  5.2.  Note  that  the  orientation  of  the  shock  in  the  Y-Z 
mesh  is  counter  to  that  of  the  numerical  domain  of  dependence  for  the 
regular  MacCormack  scheme.  By  changing  the  orientation  of  the  numerical 
domain  of  dependence  to  be  more  consistent  with  the  expected  shock 
orientation,  a  smaller  portion  of  the  numerical  domain  of  dependence  of 
a  point  near  the  shock  will  lie  across  the  discontinuity  from  that  point 
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Figure  5.2.  Typical  Alignment  of  Embedded 
Shock  in  Y-Z  Coordinate  Mesh 

The  numerical  domain  of  dependence  is  reoriented  by  reversing 
the  directions  of  the  Z-derivative  approximations  in  the  predictor  and 
corrector  stages  of  the  MacCormack  scheme;  i.e.,  backward  differences 
in  the  predictor  stage  and  forward  differences  in  the  corrector  stage. 

The  resulting  numerical  domain  of  dependence  for  this  variation  on  the 
MacCormack  scheme  is  shown  in  Figure  5.3. 

While  this  modification  to  the  numerical  scheme  cannot  eliminate 
the  oscillations  inherent  in  conservation  calculations  of  discontinuities, 
it  will,  in  some  circumstances,  decrease  the  magnitude  of  the  oscillations. 
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Figure  5.3.  Numerical  Domain  of 

Dependence  for  MacCormack 
Scheme  Inverted  in  L  Direction 

5.7  X-DIFFERENCING  FIELD  POINT  TREATMENT 

The  second  method  examined  in  this  effort  for  the  automatic 
calculation  of  embedded  shocks  is  the  X-differencing  scheme,  described 
in  Section  4.2.  In  this  approach  to  the  embedded  s  .ock  problem,  the 
governing  equations  are  written  in  non-conservation  form  and  certain 
terms  are  written  in  "split"  forms  that  allow  accurate  modeling  of  the 
domain  of  dependence.  For  example,  the  term  AP^  would  be  approximated 
as  l/2(Xz-|Pz-j  +  ^Z2^Z2^  w-*th  this  scheme,  where  the  particular  finite 
difference  expressions  used  for  Pz-j  and  Pz£  would  be  dependent  on  the 
signs  of  Xz-j  and  Xz2»  respectively. 


In  each  "split"  term,  both  the  (  )1  and  (  )2  derivatives 

are  approximated  using  a  second-order  accurate,  one-sided,  predictor- 
corrector  finite  difference  scheme.  For  either  derivative,  if  the 
corresponding  coefficient  X ^  Is  positive,  backward  differences  are 
employed;  conversely,  If  the  coefficient  is  negative,  forward  differences 
are  used.  Backward  differences  are  formed  using  Equations  (5.13)  and 
(5.11)  in  the  predictor  and  corrector  stages,  respectively,  while  forward 
differences  are  formed  using  Equations  (5.9)  and  (5.14),  as  described 
in  Section  5.2. 

It  will  be  noted  in  the  X-differencing  forms  of  the 
governing  equations,  described  in  Sections  4.2  (axi symmetric)  and  4.4 
(three  dimensional),  that  some  terms  are  not  split.  These  unsplit  terms 
are  those  that  arise  from  the  computational  transformation  to  (X,Y,Z,T) 
space  and  are  generally  small.  The  derivatives  appearing  in  these  terms 
are  approximated  using  convective  differences,  as  described  in  Section 
5.2,  maintaining  the  overall  second-order  accuracy  of  the  scheme. 


SECTION  6 


VALIDATION  OF  SOLUTION 

6.1  LIMITATIONS  Of  TECHNIQUE 

In  chis  analysis  a  time-deoendent  algorithm  has  been  de¬ 
veloped  for  the  computation  of  steady  inviscid  flows  over  ablated  re¬ 
entry  vehicle  nosetips  with  uniform  supersonic  or  hypersonic  freestream 
conditions.  This  procedure  has  been  formulated  in  a  new  coordinate 
system  that  is  capable  of  being  closely  aligned  with  any  nosetip  geometry, 
and  includes  a  technique  for  the  approximate  calculation  of  embedded 
shocks.  The  primary  limitation  that  exists  for  this  technique  is  the 
assumption  of  isentropic  flow  downstream  of  the  Dow  shock. 

Because  of  the  inviscid  flow  assumption,  the  validity  of 
this  analysis  is  limited  to  those  high  Reynolds  number  cases  where  the 
thin  boundary  layer  assumption  and  weak  interaction  theory  apply. 
Fortunately,  for  the  flight  conditions  of  reentry  vehicles  with  ablated 
nosetios,  fcr  which  this  analysis  was  undertaken,  viscous  effects  are 
generally  confined  to  a  thin  boundary  layer  adjacent  to  the  vehicle 
surface,  and  the  inviscid  flow  field  may  be  determined  independently 
of  the  boundary  layer. 

This  inviscid  assumption  fails,  however,  when  the  nosetip 
geometry  produces  separation  of  the  shock  layer  flow.  Prediction  of  the 
separated  flow  region  must  include  viscous  effects,  for  which  the  present 
analysis  is  unsuitable. 

* 

As  discussed  in  Section  1.0,  inviscid  theory  has  been 
shown  to  produce  accurate  aerodynamic  predictions  for  the  flight  conditions 
of  interest  for  all  aerodynamic  coefficients,  with  the  exception  of  the 


axial  force  coefficient,  which  can  be  significantly  affected  by  viscous 
shear  and  induced  pressure  effects.  Thus,  an  inviscid  aerodynamic  pre¬ 
diction  procedure,  such  as  developed  here,  can  be  a  valuable  tool  for 
both  the  pre-flight  and  post-flight  evaluation  efforts  relating  to  reentry 
vehicle  design  and  performance. 

Other  limitations  apply  to  the  numerical  procedure  developed 
in  this  research  that  are  conrion  to  all  numerical  fluid  flow  computa¬ 
tional  procedures.  Fundamentally,  it  is  required  for  accurate  results 
that  the  discrete  grid  points  used  in  the  numerical  calculations  be 
spaced  so  as  to  be  able  to  resolve  all  pertinent  features  of  the  body 
geometry  and  the  surrounding  flow  field.  In  particular,  a  finer  mesh 
will  be  required  in  regions  of  large  flow  gradients  to  avoid  wiggles  in 
the  numerical  solution.  (The  appearance  of  wiggles  in  ?.  numerical  solu¬ 
tion  is  indicative  of  inadequate  mesh  resolution  for  the  case  being 
computed,  or  of  some  other  error  in  the  formulation  or  application  of 
the  numerical  technique.  This  point  is  discussed  more  fully  by  Moretti32.) 
Criteria  for  the  selection  of  appropriate  mesh  spacings  for  an  earlier 
transonic  code  were  developed  by  H-'il*  Kyriss,  Truncellito,  and 
Martel! ucci 4  ;  these  criteria  are  equally  valid  for  the  procedure  developed 
in  this  report. 

Inconsistencies  between  the  mesh  point  spacing  and  flow 
gradients  can,  of  course,  be  eliminated  simply  by  using  more  mesh  points 
in  the  finite  difference  calculation.  More  core  storage  will  also  be 
required  for  the  computer  code  with  an  increased  number  of  mesh  points; 
the  only  limitations  on  the  use  of  this  analysis  in  this  regard  is  the 
available  core  storage  on  the  computer  and  the  economics  of  a  calculation 
with  a  large  number  of  mesh  points. 


6.2  CONVERGENCE  PROPERTIES 

The  time-dependent  approach  to  solving  the  steady  transonic 
problem  is  an  iterative  technique  in  which  advancement  of  the  solution 
at  any  given  iteration  is  defined  using  the  time-dependent  equations  of 
motion.  In  theory,  as  time  increases  to  infinity,  the  steady-state 
flow  field  solution  is  approached.  The  time  variable  used  In  the  time- 
dependent  technique  is  directly  analogous  to  the  iteration  number  in 
other  iterative  numerical  techniques.  In  any  iterative  technique,  it 
is  necessary  to  define  criteria  for  accepting  a  solution  as  being  con¬ 
verged.  This  section  defines  convergence  criteria  that  have  been 
developed  for  time-dependent  techniques  in  Reference  4  .  These  criteria 
are  defined  in  terms  of  fluid  dynamic  phenomena  that  arise  in  the  course 
of  calculating  the  shock  layer  about  a  body  in  supersonic  flow. 

The  criteria  detailed  below  are  sufficient  to  determine 
convergence;  i.i.t  when  they  are  satisfied,  the  convergence  of  that 
particular  solution  is  ensured.  Cases  frequently  arise,  however,  when 
a  satisfactory  solution  (from  the  standpoint  of  accurate  aerodynamic 
predictions)  can  be  obtained  without  satisfying  all  of  the  convergence 
criteria.  Acceptance  of  such  solutions  requires  judgement  on  the  part 
of  the  user.  (In  other  words,  the  criteria  presented  here  represent 
sufficient,  but  not  necessary,  conditions  for  an  acceptable  solution.) 

For  a  numerical  method  to  be  an  accurate  solution  to  a 
problem,  it  must  not  only  converge,  but  the  numerical  scheme  must  be 
consistent  with  the  problem  being  solved.  The  consistency  of  the  time- 
dependent  technique  used  in  this  analysis  is  based  on  its  discretization 
of  the  inviscid  Euler  equations,  and  its  use  of  non-dissipative 


differencing  schemes.  Thus,  the  validity  of  a  solution  obtained  with 
this  procedure  may  be  examined  solely  upon  the  degree  of  convergence 
of  the  solution  (assuming  that  the  mesh  spacing  selected  by  the  user  is 
adequate  to  define  the  details  of  the  body  geometry). 

The  convergence  criteria  developed  in  the  previous  study 
will  be  summarized  below,  and  then  discussed  in  more  detail.  These 
criteria  are: 

1.  The  "stagnation"  pressure  must  have  converged 

to  an  essentially  constant  value.  If  the  actual 
stagnation  point  is  known  to  lie  exactly  on  a 
mesh  point  in  the  finite-difference  grid,  the 
computed  value  of  stagnation  pressure  should  be 
within  0.5%  of  the  known  theoretical  value  of  pQ. 

2.  The  "stand-off  distance"  of  the  shock,  A. ,  must 

o' 

have  converged  to  a  constant  value. 

3.  The  root-mean-square  (rms)  of  the  shock  velocities, 

(ct/G)  ,  must  be  converging  (decreasing  in 
1  rms 

value)  or  have  converged  and,  in  magnitude,  must 
satisfy  the  relation  (cT/G)  s  0.004  V*,. 

4.  The  total  enthalpy  at  every  point  in  the  flow 
being  computed  must  be  within  5%  of  the  known 
steady-state  total  enthalpy. 

In  a  three-dimensional  calculation,  the  flow  stagnation 
point  does  not  usually  correspond  exactly  to  a  computational  grid  point; 
thus,  for  the  first  two  criteria  listed  above,  "p0"  is  taken  as  the 
pressure  occurring  at  the  centerline  body  point  and  "A0"  is  taken  as 
the  distance  between  the  body  point  and  the  dow  shock  point  at 


the  centerline.  To  illustrate  the  convergence  of  these  quantities, 
Figures  6.1  and  6.2  provide  examples  of  convergent  and  non-convergent 
time  histories  of  pQ. 

In  the  time-dependent  technique  on  initial  shock  shape  is 

assumed  and  is  allowed  to  adjust  its  position  during  the  coursa  of  the 

calculation.  Theoretically,  for  a  solution  that  has  converged  to  a 

steady  state,  the  shock  velocity  will  have  vanished  at  all  shock  points 

to  within  some  "epsilon"  of  the  frees tream  velocity.  On  a  more  practical 

basis,  the  criterion  described  above  ( ( cT/G)  s  0.004  V  )  is  sufficient 

1  rms 

to  determine  a  satisfactorily  converged  solution,  provided  that  the 

magnitude  of  (cT/G)  is  decreasing  from  step  to  step  (i.e.,  is  not 

diverging)when  the  criterion  is  only  marginally  satisfied.  Samples  of 

converging  and  diverging  time-histories  of  (cT/G)  are  illustrated 

1  rms 

in  the  plots  shown  in  Figures  6.3  and  6.4.  One  caveat  is  required 

in  the  assessment  of  (ct/G)  .  Since  (cT/G)  is,  in  a  sense,  an 

'  rms  1  rms 

average  of  individual  shock  velocities,  it  is  possible  that  the  criterion 

on  (ct/G)  be  satisfied,  while  one  or  two  individual  shock  velocities 
1  rms 

are  relatively  large.  Thus,  judicious  assessment  of  convergence  based 

on  shock  velocity  should  include  examination  of  the  individual  shock 

velocities  as  well  as  the  value  of  (cT/G) 

1  ms 

The  final  criterion  for  convergence  is  based  upon  the 
conservation  of  total  enthalpy.  Since  the  total  enthalpy  is  not  con¬ 
stant  in  an  unsteady  flow,  an  indication  of  the  convergence  and  accuracy 
of  a  time-deper.dent  solution  can  be  obtained  by  examining  the  difference 
between  the  computed  total  enthalpy  at  each  point  and  the  known  steady 
state  value  (which  is  equal  to  the  freestream  total  enthalpy). 
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Figure  6.2.  Typical  p0  History  (Not  Converged) 
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Figure  6,3.  Typical  (ct/g)  History  (Converged) 
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Defining 


AH  =  (H-HJ/H^  (6.1) 

a  suitable  criterion  for  convergence  has  been  found  to  be 

| AH |  s  0.05  .  (6.2) 

The  question  of  convergence  of  the  time-depende.v;  solution 
can  be  sunmarized  as  follows:  when  a  time-dependent  solution  satisfies 
the  convergence  criteria  described  above,  it  can  be  used  with  confidence 
to  provide  accurate  and  reliable  aerodynamic  predictions,  provided  that 
the  mesh  adequately  resolves  the  important  details  of  the  body  geometry 
and  of  the  surrounding  flow  field. 

6.3  VALIDATION  OF  NOSETIP  SOLUTION 

In  this  section  the  numerical  procedures  developed  in  this 
research  are  validated  by  comparisons  of  predictions  to  experimental  data 
and,  where  appropriate,  predictions  obtained  with  other  numerical  tech¬ 
niques.  The  validation  process  described  in  this  section  pertains  only 
to  the  nosetip  solution  procedures;  the  ability  of  these  transonic 
procedures,  when  coupled  to  a  steady  supersonic  afterbody  code,  to  make 
accurate  determinations  of  total  vehicle  inviscid  aerodynamics  is 
demonstrated  in  Section  6.4. 

Two  versions  of  the  three-dimensional  time-dependent 
inviscid  code  formulated  in  the  new  coordinate  system  (based  on  a  series 
of  conformal  transformations)  have  been  developed  in  this  effort.  The 
basic  version,  which  is  denoted  by  CM3DT(NC),  is  formulated  in  terms  of 
the  non-conservation  dependent  variables  and  is  incapable  of  treating 


123. 


embedded  shocks.  The  version  based  on  the  X-differencing  scheme, 
providing  a  method  for  treating  embedded  shocks,  is  denoted  by  CM3DT(X). 

The  first  step  in  the  validation  process  is  the  demonstra¬ 
tion  of  the  ability  of  the  new  coordinate  system  to  permit  accurate 
calculations  of  the  inviscid  flow  over  a  wide  variety  of  body  shapes, 
including  shapes  that  could  not  previously  be  treated  with  transonic 
codes  formulated  in  standard  coordinate  systems  (e.g.,  spherical).  Use 
of  this  new  coordinate  system  does  not,  however,  reduce  the  accuracy  of 
calculations  for  shapes  that  are  well-suited  to  analysis  using  a  standard 
coordinate  system. 

To  illustrate  the  ability  of  the  new  coordinate  system  to 
treat  shapes  that  are  well-suited  to  calculations  using  a  standard 
coordinate  system,  Figures  6.5  and  6.6  present  predictions  of  the  bow 
shock  shape  and  surface  pressure  distribution  for  a  sphere  in  equilibrium 
air  at  an  altitude  of  100  KFT  with  a  freestream  velocity  of  20,000  ft/sec. 
Two  predictions  are  made  with  the  CM3DT(NC)  technique,  using  both  a 
coarse  mesh  (6  x  9;  i.e.,  six  points  across  the  shock  layer  and  nine 
points  along  the  body)  and  a  fine  mesh  (11  x  17).  The  other  predictions 
shown  for  this  flow  have  been  obtained  from  the  inverse  technique  of 
Lomax  and  Inouye33  and  from  the  technique  developed  by  Kyriss  ?nd 
Harris0  ,  which  is  an  explicit  time-dependent  finite-difference  code 
formulated  in  a  spherical  coordinate  system.  (The  method  of  Kyriss  and 
Harris  has  been  extensively  validated  by  comparisons  to  experimental 
data,  e.g.,  by  Hall,  Kyriss,  Truncellito,  and  Martel! ucci  4  and  by 
Hall  and  Nowlan5  .) 
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Figure  6,5,  Shock  Shape  Predictions  for  a 
Sphere  in  Equilibrium  Real  Gas 
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The  results  obtained  with  the  CM3DT(NC)  code  show 
excellent  agreement  with  the  other  predictions  for  th's  case  and, 
furthermore,  are  also  seen  to  be  invariant  with  refinement  of  the  mesh, 
verifying  the  consistency  of  the  numerical  approximations.  This  case 
also  serves  to  illustrate  the  real  gas  thermodynamic  capabilities  of 
the  CM3DT  procedure. 

The  ability  of  the  new  coordinate  system  to  treat  less 
regular  shapes  is  shown  in  Figures  6.7-6.10,  which  present  comparisons 
of  bow  shock  shape  and  surface  pressure  distribution  predictions  for 
the  axi symmetric  PANT  Triconic  (described  by  Jackson  and  Baker  in 
Reference  21)  to  experimental  data  at  =  5  and  a  B  0°.  (Although 
an  indented  shape,  no  embedded  shock  forms  on  this  configuration  at 
these  flow  conditions.)  Figure  6.7  presents  the  surface  pressure  distri¬ 
bution  predicted  with  the  CM3DT(NC)  procedure  for  this  configuration, 
which  has  a  small  radius  corner  at  the  shoulder  leading  back  to  an  aft 
cone,  as  evident  in  Figure  6.10,  The  agreement  between  the  prediction 
and  the  data  is  seen  to  be  good  in  Figure  6.7,  except  for  oscillations 
arising  in  the  vicinity  of  the  corner.  These  oscillations  are  the  result 
of  inadequate  resolution  of  the  finite  difference  grid  in  the  vicinity 
of  the  corner,  as  discussed  and  illustrated  by  Hall,  Kyriss,  Truncellito, 
and  Martellucci  l;  ,  The  influence  of  this  sharp  corner  is  made  more 
evident  by'  the  results  obtained  from  a  similar  calculation,  in  which 
the  sharp  corner  and  cone  were  removed.  As  shown  in  Figure  6.8,  no 
oscillations  appear  in  the  CM3DT(NC)  prediction  of  surface  pressure  when 
the  sharp  corner  is  eliminated. 
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Figure  6,7.  CM3DT  (NC)  Prediction  of  Surface 
Pressure  Distribution  for  the 
PANT  Triconic  (with  cone)  at 
fL  =  5,  a  =  0° 
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Figure  6.8.  CM3DT  (NC)  Prediction  of  Surface 
Pressure  Distribution  for  the  PANT 
Trtconic  (without  cone)  at  fl,  =  Sj 

a  =  0° 
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Pressure  Distribution  for  the  PANT 
Triconic  (with  cone)  at  M.  =  5,  <»* 


Figure  6.10,  CM3DT  Shock  Shape  Prediction  for  yhe  PANT  Triconic  at  ft  = 


It  was  discovered  during  the  course  of  this  research  that 
the  X-differencing  scheme,  in  addition  to  its  shock-capturing  capability, 
has  the  additional  ability  of  eliminating  the  spurious  oscillations  that 
arise  at  such  sharp  corners.  This  capability  is  illustrated  in  Figure 
6.9,  which  depicts  the  surface  pressure  distribution  for  the  PANT  Triconic 
(with  the  sharp  corner  and  conical  aft  section  included)  as  predicted 
by  the  CM3DT(X)  procedure.  The  lack  of  oscillations  in  this  calculation 
is  a  manifestation  of  the  ability  of  the  X-differencing  scheme  to 
accurately  model  the  physical  domain  of  dependence,  preventing  oscillations 
from  propagating  upstream  in  supersonic  flow. 

The  final  comparison  for  the  PANT  Tri conic  predictions 
is  shown  in  Figure  6.10,  demonstrating  the  agreement  between  predictions 
and  experimental  data  for  the  bow  snock  shape.  (All  CM3DT  predictions, 
both  with  and  without  the  sharp  corner,  for  the  bow  shock  shape  are 
essentially  equivalent.) 

It  should  be  noted  that  attempts  at  computing  this  slender 
shape  with  the  technique  of  Kyriss  and  Harris  6  ,  which  is  formulated  in 
a  spherical  coordinate  system,  were  unsuccessful,  because  of  the  inability 
of  the  spherical  coordinate  system  to  be  closely  aligneo  with  the  body 
geometry . 

Further  evidence  of  the  abilities  of  the  new  coordinate 
system  is  provided  in  Figure  6.11,  which  presents  a  comparison  of  surface 
pressure  predictions  and  experimental  data  for  the  PANT  Simple  Biconic 
(described  by  Jackson  and  Baker21)  at  =5  and  a  =  5°.  This  axisym- 
metric  configuration  is  a  45°  sphere-cone  nosetip,  with  a  rounded 
shoulder  leading  to  a  6°  aft  cone.  As  seen  in  this  figure,  the 
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SURFACE  PRESSURE  p/p 


Figure  6.11.  Surface  Pressure  Predictions  for  the 
PANT  Simple  Biconic  at  Mw  =  5,  <*  =  5° 


predictions  of  both  CM3DT(NC)  and  the  technique  of  Kyriss  and  Harris8 
show  good  agreement  with  the  data,  except  for  the  oscillations  created 
by  the  lack  of  adequate  resolution  at  the  shoulder.  Not  apparent  from 
this  figure,  however,  is  the  improved  convergence  behavior  obtained 
with  the  new  transonic  technique.  After  400  time  steps,  using  the 
same  number  of  grid  points,  and  computing  approximately  the  same  elapsed 
physical  time,  the  root-mean-squares  of  the  shock  velocities,  which, 
as  discussed  in  Section  6.2,  are  reliable  indicators  of  convergence, 
differ  by  an  order  of  magnitude.  For  CM3DT(NC),  (cT/G)  /V  =  0.004, 
while  the  equivalent  quantity  for  the  other  technique  is  0.056,  which 
does  not  satisfy  the  convergence  criterion. 

Another  good  indicator  of  convergence  for  a  time-dependent 
technique  is  the  maximum  error  in  computed  total  enthalpy.  For  CM3DT(NC), 
the  maximum  error  is  8.6%  for  this  case,  compared  to  a  maximum  error  of 
28.5%  for  the  technique  of  Kyriss  and  Harris,  further  illustrating  the 
benefits  of  the  new  coordinate  system. 

The  ability  of  the  CM3DT(NC)  code  to  treat  asymmetric  shapes 
at  angle  of  attack  using  the  new  coordinate  system  is  illustrated  in 
Figure  6,12,  showing  predictions  of  the  bow  shock  shape  for  the  Blunt-1 
shape  at  =  13,4  and  a  *  3°.  (This  asymmetric  nosetip  shape  was 
derived  in  a  nosetip  reconstruction  analysis  described  by  Hall  and 
Nowlan5  .)  The  predicted  shock  shape  is  seen  to  agree  well  with  that 
predicted  using  the  technique  of  Kyriss  and  Harris8  . 

The  final  step  in  the  validation  process  for  the  nosetip 
solution  is  the  demonstration  of  the  ability  of  the  Cti3DT(X)  procedure 
to  predict  flows  over  shapes  that  produce  an  embedded  shock.  In  Figure 


134. 


6,12.  Bow  Shock  Shape  Predictions  for  the 
the  Blunt -1  Shape  at  =  13.4,  a  = 


6.13,  a  comparison  is  shown  of  the  predicted  surface  pressure  distri¬ 
bution  and  the  experimental  data  for  the  PANT  Triconic  at  M  =5  and 

00 

a  *  10°;  the  agreement  is  seen  to  be  excellent.  Schlieren  photographs 
of  this  configuration  at  this  angle  of  attack,  which  may  be  found  in 
Reference  21,  clearly  indicate  the  presence  of  an  embedded  shock  in 
the  lee  plane.  This  comparison  serves  to  illustrate  the  capability  of 
the  three  dimensional  version  of  the  A-differencing  scheme  to  compute 
inviscid  nosetip  flow  fields  with  embedded  shocks.  (The  capabilities  of 
the  axisymmetric  version  of  the  A-differencing  scheme  have  been  demon¬ 
strated  in  Section  4.) 

6.4  PREDICTION  OF  TOTAL  VEHICLE  INVISCID  AERODYNAMICS 

In  the  preceding  section  the  ability  of  the  CM3DT  technique  to 
accurately  predict  inviscid  flow  fields  over  ablated  reentry  vehicle 
nosetips  has  been  demonstrated.  To  complete  the  validation  process  for 
this  technique,  it  remains  only  to  demonstrate  the  capability  of  the 
CM3DT  nosetip  code,  when  coupled  with  an  existing  supersonic  afterbody 
code,  to  make  accurate  predictions  of  total  vehicle  inviscid  aerodynamics. 
For  this  validation  effort,  the  CM3DT  code  has  been  coupled  to  the  super¬ 
sonic  afterbody  code  of  Kyriss  and  Harris8  ,  which  is  a  steady,  forward 
marching  finite  difference  technique  formulated  in  a  cylindrical  co¬ 
ordinate  system.  (Coupling  of  the  nosetip  and  afterbody  codes  requires 
interpolation  on  the  nosetip  solution  to  determine  the  flow  variables 
in  the  initial  data  plane  of  the  afterbody  solution  and  an  integration 
to  determine  the  forces  and  moments  acting  on  the  nosetip.) 
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Figure  6,13.  Surface  Pressure  Predictions  for  the 
PANT  Triconic  at  =  5,  «  =  10° 
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In  the  calculation  of  total  vehicle  aerodynamics  it  is  of 
critical  importance  to  accurately  model  the  influence  of  the  nosetip 
shape  on  the  afterbody  flow  field,  and  thus  on  the  resulting  forces 
and  moments  acting  on  the  vehicle.  For  example,  the  importance  of 
the  downstream  influence  of  an  asymmetric  nosetip  in  the  determination 
of  the  trim  angle  of  attack  (cty)  for  a  ballistic  reentry  vehicle  has 
been  demonstrated  by  Hall  and  Nowlan5  .  Direct  coupling  of  numerical 
flow  field  calculations  for  nosetips  with  afterbody  procedures  will 
automatically  include  the  nosetip's  influence  on  the  afterbody  flow 
field. 

The  first  validation  case  for  the  CM3DT  code  coupled  to  the 
afterbody  solution  is  presented  in  Figure  6.14.  This  figure  depicts  a 
comparison  of  predictions  of  the  normal  force  coefficient  (CN)  and  the 
pitching  moment  coefficient  (Cm)  as  a  function  of  angle  of  attack  for 
a  9°  cone  with  a  spherical  nose,  with  a  vehicle  bluntness  ratio  (RN/Rg) 
of  15%,  at  =  20  in  ideal  gas  (y  =  1.4)  with  no  sideslip  (3  =  0°). 

The  prediction  labeled  CM3DT  was  obtained  using  the  CM3DT(NC)  code  to 
provide  the  initial  data  required  for  the  Kyriss  and  Harris8  after¬ 
body  code;  the  other  prediction  was  obtained  using  an  axisymmetric 
calculation  of  the  spherical  nosetip  flow  field  (in  wind-fixed  coordi¬ 
nates),  with  the  required  initial  data  being  obtained  by  suitable  rotations 
of  the  spherical  solution,  as  described  in  Reference  8  .  The  accuracy 
of  the  Kyriss  and  Harris  afterbody  code,  and  of  the  Kyriss  and  Harris 
axisymmetric  and  three-dimensional  transonic  codes  for  blunt,  convex 
nosetips,  has  been  extensively  demonstrated  through  comparisons  of 
predictions  to  wind  tunnel  data,  as  shown,  for  example,  in  References 
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Figure  6.14.  Normal  Force  and  Pitching  Moment 

Coefficients  vs.  Angle  of  Attack  for 
9°  Sphere-Cone  of  15%  Bluntness  at 
H.  -  20,  0  =  0° 


4,5,  and  8  .  As  is  evident  in  Figure  6.14,  the  agreement  between 
the  predictions  of  C^(a)  and  Cm(a)  is  excellent. 

The  comparison  of  aerodynamic  predictions  for  spherically 
blunted  cones  obtained  with  these  two  techniques  is  continued  in  Figure 
6.15.  This  figure  illustrates  the  excellent  agreement  obtained  between 
the  two  predictions  of  pitch  center  of  pressure  as  a  function  of 
bluntness  ratio  for  the  9°  sphere-cone  at  *  20  and  a  =  5°.  (The 
pitch  center  of  pressure  is  defined  as  xcp/L ^  =  xcg/L^  -  Cm/CN  when 
the  reference  length  used  in  the  non-dimensional ization  of  the  pitching 
moment  is  the  virtual  cone  length.) 

The  ability  of  CM3DT  to  provide  accurate  initial  data  to 
the  afterbody  calculation  for  spherically-nosed  vehicles  is  further 
demonstrated  in  Figure  6.16.  This  comparison  is  similar  to  that  shown 
in  Figure  6.14,  in  that  predictions  are  obtained  for  CN(a)  and  Cm(a) 
for  a  15%  blunt  9°  sphere-cone  at  =  20,  except  that  the  vehicle  is 
at  k  onstant  sideslip  angle  (B)  of  5°.  The  agreement  between  the  two 
predictions  is  again  seen  to  be  excellent,  and  serves  to  verify  the 
ability  of  the  CM3DT  to  compute  flow  fields  on  nosetips  in  both  pitch 
and  yaw. 

Figures  6.17  and  6.18  present  a  comparison  of  inviscid  aerodynamic 
.‘on?  ained  for  a  vehicle  with  a  blunt  asymmetric  nosetip.  The 
vehicle  is  a  6°  cone  with  a  nominal  bluntness  ratio  of  25%;  the  nosetip 
is  the  Blunt-1  shape  (illustrated  in  Figure  6.12),  which  was  derived  in 
a  nosetip  shv  . econstruction  effort  described  by  Hall  and  Nowlans  . 

(In  Reference  5  ,  the  Blunt-1  shape  was  selected  as  a  plausible  nosetip 
shape  for  an  actual  flight  vehicle  at  20  KFT,  based  on  actual  recession 
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Figure  6.15.  Pitch  Center  of  Pressure  vs.  Bluntness 
Ratio  for  9°  Sphere-Cone  at  PL  =  20, 
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Figure  6.16.  Normal  Force  and  Pitching  Moment 
Coefficients  vs.  Angle  of  Attack 
for  9°  Sphere-Cone  of  15%  Bluntness 
at  =  20,  b  *  5° 
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Figure  6.17.  Normal  Force  Coefficient  vs.  Angle 
of  Attack  for  BLUNT-1  Configuration 
at  M.  =  13,4,  20  KFT  Altitude 
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Figure  6.13,  Pitching  Moment  Coefficient  vs,  Angle 
of  Attack  for  BLUNT-1  Configuration 
at  =  13.4 ,  20  KFT  Altitude 


measurements  and  on  the  consistency  of  aerodynamic  characteristics  derived 
from  motion  data  with  those  computed  with  the  transonic  nosetip  and 
supersonic  afterbody  codes  of  Kyriss  and  Harris8  .  The  computed  trim 
angle  of  attack,  *  3.35°,  was  consistent  with  that  measured  by  an 
onboard  magnetometer,  as  indicated  in  Figure  6.18.) 

A  comparison  of  C^(a)  predictions  for  this  configuration  is 
presented  in  Figure  6.17,  and  fair  agreement  is  evident  between  the 
predictions  obtained  with  the  CM3DT(NC)  technique  and  with  the  three- 
dimensional  transonic  code  of  Kyriss  and  Harris.  Similar  agreement  is 
seen  in  the  Cm(a)  predictions  shown  in  Figure  6.18.  The  CM3DT  pre¬ 
dictions  indicate  a  trim  angle  of  attack  of  3.58°,  which  is  within  the 
Oj  range  measured  by  the  magnetometer. 

The  above  comparisons  have  demonstrated  the  ability  of  the 
CM3DT  code  to  produce  nosetip  flow  field  solutions  and  afterbody  solution 
initial  data  for  spherical  and  blunt  convex  nosetip  shapes  with  an  accuracy 
equivalent  to  that  of  the  extensively  validated  crar.sonic  technique  of 
Kyriss  and  Harris8  .  Because  the  CM3DT  code  has  been  developed  to 
extend  the  range  of  nosetip  geometries  for  which  accurate  numerical 
inviscid  aerodynamic  predictions  can  be  made,  however,  the  remaining 
validation  cases  to  be  documented  in  this  section  must  rely  on  compari¬ 
sons  of  CM3DT  solutions  with  wind  tunnel  data  for  complete  reentry 
vehicles,  rather  than  with  other  numerical  solutions. 

The  first  case  to  be  considered  in  this  comparison  of  numerical 
predictions  and  experimental  data  is  the  N8  nosetip,  shown  in  Figure 
6.19,  mounted  on  a  6.3°  cone,  with  a  nominal  bluntness  of  25%. 
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Experimental  data  were  obtained  at  Mach  8  in  Tunnel  B  at  the  Arnold 
Engineering  Development  Center  (AEDC);  details  of  the  tests  and  the 
resulting  data  are  described  by  Hahn  and  Little34.  The  N8  nosetip  is 
a  "tilted  cone"  shape,  with  relatively  small  shoulder  radii,  and  is 
intended  to  simulate  an  asymmetric  shape  that  could  result  from 
turbulent  ablation  on  the  nose.  Because  of  the  fairly  sharp  corners 
on  this  shape,  the  CM3DT(X)  code  was  used  to  generate  the  nosetip  flow 
field  solution,  since  the  X-differencing  scheme  has  been  found  to 
eliminate  wiggles  arising  at  sharp  corners,  as  discussed  in  Section  6.3. 

Comparisons  of  the  predictions  of  the  normal  force  and  pitching 
moment  coefficients  to  the  experimental  data  are  shown  in  Figures  6.20 
and  6.21,  respectively.  These  comparisons  show  good  agreement  between 
the  predictions  and  the  data.  Also,  the  CM3DT  technique  has  accurately 
computed  the  trim  angle  of  attack  to  be  3.86°  (based  on  a  linear  inter¬ 
polation  between  values  of  Cm  computed  at  a  =  2°  and  a  =  4°),  compared 

to  the  experimentally  determined  value  of  3.52°  (based  on  a  linear 
interpolation  between- data  points  at  a  =  3°  and  a  =  4°). 

The  final  validation  cases  to  be  presented  in  this  section 
involve  comparisons  between  predictions  and  tunnel  measurements  of 
aerodynamic  forces  and  moments  for  a  6°  cone  (25%  nominal  bluntness) 
with  two  axisymmetric  indented  nosetips.  The  tests  were  conducted  in 

Tunnel  F  at  AEDC  at  a  Mach  number  of  11.6  and  are  described  by  Boudreau, 

Crain,  and  Edenfield3.5  The  nosetips  tested  were  the  Very  Mildly  In¬ 
dented  Body  (VMIB)  and  the  Mildly  Indented  Body  (MIB),  described  by 
Reeves,  Todisco,  Lin,  and  Pallone20,  and  illustrated  in  Figures  4.4 
and  4.7,  respectively. 
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Since  both  the  VMIB  and  MIB  nosetips  are  indented  and  test"  were 

conducted  at  non-zero  angles  of  attack,  it  was  necessary  to  use  the 
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CM3DT(A)  code  because  of  the  shock-capturing  capability  of  the  A-dif- 
ferencing  scheme  for  the  treatment  of  embedded  shocks.  Figures  6.22  and 
6.23  present  comparisons  of  C^(a)  and  Cm(a)  predictions  to  the  data, 
respectively,  for  the  VMIB  nosetip,  and  the  agreement  is  seen  to  be 
excellent. 

The  same  comparisons  are  presented  in  Figures  6.24  and  6.25  for 
the  vehicle  with  the  MIB  nosetip.  Again,  the  overall  agreement  between 
the  predictions  and  the  data  is  seen  to  be  good.  (Discrepancies  are 
apparent,  however,  between  the  predicted  and  measured  values  of  and  Cm 
at  a  =  4°  for  this  configuration.  The  validity  of  the  CM3DT  calculation 
at  this  angle  of  attack,  as  measured  by  the  convergence  of  the  computation, 
is  comparable  to  that  of  the  calculations  at  other  angles  of  attack,  where 
good  agreement  between  predictions  and  data  is  evident.) 

Further  validation  of  the  VMIB  and  MIB  calculations  is  presented 
in  Figure  6.26,  which  depicts  the  predicted  pitch  center  of  pressure  loca¬ 
tions  compared  to  the  experimental  center  of  pressure  data  presented  in 
Reference  8.  The  agreement  is  seen  to  be  excellent. 

In  conclusion,  the  ability  of  the  CM3DT(A)  code  to  compute  in- 
viscid  flow  fields  on  indented  nosetips  at  angle  of  attack,  when  coupled 
with  a  supersonic  afterbody  code,  has  resulted  in  a  technique  that  is 
capable  of  producing  accurate  aerodynamic  predictions  for  reentry  vehicles 
with  realistic  ablated  nosetip  shapes,  as  demonstrated  above.  This  capa¬ 
bility  represents  a  significant  extension  of  the  applicability  of  numerical 
techniques  to  the  evaluation  of  reentry  vehicle  inviscid  aerodynamics. 
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NORMAL  FORCE  COEF 


ANGLE  OF  ATTACK  (DEG.) 

Figure  6.22,  Normal  Force  Coefficient  vs.  Angle  of 

Attack  6or  VMIB  Configuration  at  =  11.6 
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O  EXPERIMENTAL  DATA^ 
Re^  =  9  x  106/ft 

-  CM3DT  (X) 


ANGLE  OF  ATTACK  (.DEG. ) 


Figure  6,24,  Normal  Force  Coefficient  vs.  Angle  of 

Attack  for  MIB  Configuration  at  PL  *  11.6 


Figure  6,25.  Pitching  Moment  Coefficient  vs,  Angle 
of  Attack  for  MIB  Configuration  at  M« 
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SECTION  7 


CONCLUSIONS 

The  goal  of  this  research  effort  has  been  the  development  of 
a  computational  technique  for  the  prediction  of  inviscid  flow  fields 
about  ablated  reentry  vehicle  nosetips.  The  development  of  this  capa¬ 
bility  has  required  the  elimination  of  two  major  deficiencies  in 
existing  supersonic  blunt  body  flow  field  procedures.  First,  a  co¬ 
ordinate  system  hus  been  developed  that  is  capable  of  being  closely 
aligned  with  ablated  nosetip  shapes,  producing  coordinate  surfaces 
that  are  either  nearly  parallel  or  nearly  normal  to  the  nosetip  surface. 
Second,  an  approximate  capability  has  been  developed  for  the  calcula¬ 
tion  of  embedded  shocks  on  indented  nosetips. 

The  new  three-dimensional  coordinate  system  is  based  on 
sequences  of  conformal  transformations  that  are  carried  out  independently 
in  each  meridional  plane.  The  conformal  transformations  developed  for 
this  effort  are  defined  in  terms  of  "hinge  points",  which  are  discrete 
points  selected  such  that  the  body  contours  are  modeled  in  an  approxi¬ 
mate  manner. 

A  three-dimensional  blunt  body  code  has  been  developed  using 
this  new  coordinate  system,  with  the  governing  equations  written  in 
non-conservation  form.  A  second-order  accurate,  explicit  procedure 
has  been  used  to  integrate  the  time-dependent  equations,  with  the 
steady  state  solution  being  sought  as  the  asymptotic  limit  of  an 
unsteady  flow.  Body  points  and  bow  shock  points  are  computed  using 
special  algorithms  that  are  based  on  the  method  of  characteristics. 
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This  new  three-dimensional  blunt  body  code  has  extended  the 
range  of  nosetip  geometries  for  which  inviscid  flows  can  be  success¬ 
fully  computed.  For  ablated  nosetlps  without  embedded  shocks,  this 
technique  has  demonstrated  excellent  agreement  with  experimental  data 
and  (for  sufficiently  regular  shapes)  with  other  numerical  techniques. 

As  a  first  step  in  treating  the  embedded  shock  problem  on 
ablated  noset ips,  two  shock-capturing  techniques  have  been  examined 
for  the  axi symmetric  problem.  These  two  procedures,  the  conservation 
formulation  and  the  X-differencing  approach,  can  only  approximate 
discrete  embedded  shocks,  but  do  not  require  special  logic  to  detect 
the  embedded  shocks  or  to  treat  their  movement  through  the  coordinate 
mesh  during  the  time-dependent  calculation.  Comparisons  of  these  two 
techniques  have  been  carried  out  by  developing  two  axisynmetric  time- 
dependent  blunt  body  codes,  using  the  new  coordinate  system  developed 
in  the  first  part  of  this  effort. 

The  conservation  formulation,  in  which  the  dependent  variables 
used  are  derived  from  the  integral  conservation  laws,  was  found  to 
require  numerical  damping  to  control  the  oscillations  produced  when 
capturing  a  shock.  This  artificial  damping  has  been  demonstrated  to 
have  the  potential  of  distorting  the  shock  layer  being  computed, 
leading  to  smooth  but  inaccurate  results. 

In  the  X-differencing  approach,  in  which  the  non-conservation 
form  of  the  governing  equations  is  solved  using  a  finite  difference 
scheme  that  accurately  models  the  domain  of  dependence  of  each  mesh 
point,  no  numerical  damping  was  required.  In  general,  the  X-differencing 
scheme  has  been  found  to  produce  better  agreement  with  experimental  data 
for  nosetips  with  embedded  shocks  than  the  conservation  approach. 
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The  X -differencing  scheme  is  at  best,  however,  an  approximate  embedded 
shock  solution,  since  there  is  no  mechanism  for  an  increase  in  entropy 
across  the  "captured"  embedded  shock. 

Based  on  this  comparison  of  axi  syrnnetri  c  shock-capturing 
techniques,  the  X-differencing  scheme  was  selected  for  extension  to 
the  three-dimensional  problem,  using  the  new  coordinate  system.  (This 
effort  is  the  first  use  of  the  X-differencing  scheme  for  a  three- 
dimensional,  time-dependent  problem.)  The  three-dimensional  X-differencing 
scheme  has  produced  good  agreement  with  experimental  data  for  indented 
nosetips  at  angle  of  attack. 

In  addition  to  its  shock-capturing  abilities,  the  X-differencing 
scheme  has  also  demonstrated  an  ability  to  eliminate  oscillations  that 
appear  in  other  blunt  body  solutions  at  sharp  corners  in  supersonic  flow. 
This  capability  arises  from  the  accurate  modeling  of  the  domain  of 
dependence  of  each  grid  point  by  the  X-differencing  scheme,  preventing 
disturbances  from  propagating  upstream  in  supersonic  flow.  The  capa¬ 
bilities  of  the  X-differencing  scheme  are  evidence  of  the  importance 
of  considering  the  physics  of  the  flow  when  developing  a  numerical 
simulation. 

The  codes  developed  in  this  effort  have  significantly  in¬ 
creased  the  range  of  ablated  nosetip  geometries  for  which  inviscid 
aerodynamic  predictions  can  be  obtained.  Coupling  these  nosetip  codes 
to  an  existing  supersonic  afterbody  code  has  provided  a  unique  capa¬ 
bility  for  the  determination  of  inviscid  aerodynamic  performance 
for  an  entire  reentry  vehicle,  both  for  pre-flight  aerodynamic 
predictions  ^.nd  for  post-flight  performance  analyses.  The  utility 
of  these  codes  is  enhanced  by  their  relative  efficiency,  with  a 
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typical  nosetip  flow  field  calculation  requiring  approximately  three 
minutes  on  the  CDC  Cyber  176  computer*  allowing  the  use  of  these 
codes  on  a  routine  basis  for  design  and  evaluation  efforts. 

Extension  of  the  current  effort  would  require  the  develop¬ 
ment  of  a  procedure  for  treating  embedded  shocks  as  discrete  dis¬ 
continuities.  Such  a  shock-fitting  technique,  in  which  the  Rankine- 
Hugoniot  conditions  are  strictly  enforced  across  the  embedded  shocks, 
would  remove  the  approximations  Inherent  in  the  isentropic  X-differencing 
scheme.  The  development  of  a  shock-fitting  algorithm  will  be  simplified 
by  the  use  of  the  X-differencing  scheme  since  the  X-scheme  provides 
a  convenient,  reliable  method  for  the  detection  of  shocks.  Coupled 
with  the  generalized  coordinate  system  developed  in  this  effort,  a 
shock-fitting  technique  would  represent  the  ultimate  capability  for 
flow  field  predictions  for  ablated  reentry  vehicle  nosetips,  subject  to 
the  limitations  inherent  in  the  assumption  of  inviscid  flow. 
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APPENDIX 


COEFFICIENTS  FOR  THE  SHOCK  ACCELERATION  EQUATION 

In  Section  5.4,  a  form  of  the  characteristic  compatibility 
condition  for  bow  shock  points  written  in  terms  of  the  shock  accelera¬ 
tion  Cyj  is  given  by  Equation  (5.36).  Because  of  the  complexity  of 
the  resulting  equation,  the  shock  acceleration  is  written  for  con¬ 
venience  in  terms  of  a  number  of  coefficients,  which  are  derived  in 
this  Appendix. 

First,  it  is  necessary  to  find  the  time  variation  of  the 
unit  shock  normal  by  differentiating  Equations  ( 5 . 29 ) - ( 5 . 31 ) ,  re¬ 
sulting  in 

N2t  =  -N2(4>2ct  +  vT/v)  (A.l) 

Niy  *  -  c^N2y  -  Ng  C^y  (A. 2) 

N3t  =  -N3(CCy/Gy  +  Vy/v)  +  [Cy(G^/G  +  c^) 

•  Vet  '  ceT,/yu  (A-3> 

where 

Vy  =  Ng  [G{-<j)2Cy  (1+C^  )  +  C^C^y} 

~^ct  “  Ce)2/G2y3 

+  ct^ G4>/G  +  "  ^<}>C5T  "  c9T^n<J)”^4»cC”c0^Gy2'*  ‘ 

Defining 

Uccy  =  C]  Cyy  +  C2  (A.  5) 
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and  differentiating  Equation  (5.24)  results  in 


c1  =  -n2/g 

^2  =  ^1CT^“V®^1  +  u»N2)  ~  cj24>2^2^ 

(A.6) 

+  uqoN1t  +  ‘  ct/g)  n2j  +  wooN3T 

(A. 7) 

where  it  has  been  noted  that 

U«T  =  - 

V~T  =  u®^icT 

Other  required  time  derivatives  are  defined  as 

PT  =  C3  CTT  +  C4 

(A. 8) 

UT  =  C5  CTT  +  C6 

(A. 9) 

which  are  to  be  evaluated  from 


p  =  3P  q 

T 

(A. 10) 

Su 

UT  "  3S»U»T 

(A. 11) 

resulting  in 

C  ■  C 

3  3u“  H 

00 

(A. 12) 

C  =  — Jc—  c 

4  lul 

00 

(A. 13) 

c  =  lu_  c 

5  3G  H 

00 

(A. 14) 

c  =  c 

L6  3T  l2 

(A. 15) 
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The  derivatives  SP/au^  and  30/30^  may  be  evaluated  by  straightforward 

differentiation  of  the  well-known  Rankine-Hugoniot  conditions  for  an 

ideal  gas;  for  equilibrium  real  gas  thermodynamics,  where  the  Rankine- 

Hugoniot  conditions  must  be  solved  by  iteration,  these  derivatives  may 

be  evaluated  numerically  from  tables  of  shock  properties  as  functions 

of  u  . 

00 

Velocity  components  downstream  of  the  shock,  given  by 
Equations  (5.34)-(5.36)  may  be  differentiated  to  give 


UT  =  ^9  CTT  +  ^10  (A. 16) 

Vj  =  C-j i  Cj-j-  +  (A.  17) 

WT  =  C 1 3  CTT  +  C14  (A. 18) 

where 

Cg  =  (C5  -  N1  (A. 19) 

C10  =  "v«<hcT  +  (c6  '  c2^  Nl 

+  (u  -  uj  N1t  (A. 20) 

^11  =  ^5  ”  ^l  ^  ^2  (A. 21) 

C12  =  u«^ict  +  (C6  "  ^  ^2 

+  (u  -  uj  N2t  (A. 22) 

C13  =  ^C5  "  Cl^  N3  (A. 23) 

C-J4  =  (Cg  -  Cg)  Ng  +  (u  -  Q.)  N^j  .  (A. 24) 
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5301  Bolsa  Avenue 
Huntington  Beach,  CA  92647 

PDA  Engineering 
Attn:  M.  Sherman 
(3)  1560  Brookhollow  Drive 
Santa  Ana,  CA  92705 

Sandia  Laboratories 
P.  0.  Box  5800 
Attn:  Library 
Albuquerque,  NM  87115 

Science  Appl ications,  Inc. 

Attn:  A.  Martel  luce i 
994  Old  Eagle  School  Road 
Suite  1018 
Wayne,  PA  19087 


